--- /dev/null
+C Program Mult_Gen
+ SUBROUTINE multgen
+ implicit none
+
+CCC Documentation and Description:
+CCC
+C This code is intended to provide a quick means of producing
+C uncorrelated simulated events for event-by-event studies,
+C detector acceptance and efficiency studies, etc. The
+C user selects the number of events, the one-particle distribution
+C model, the Geant particles to include, the ranges in transverse
+C momentum, pseudorapidity and azimuthal angle, the mean
+C multiplicity for each particle type for the event run, the
+C mean temperature, Rapidity width, etc., and the standard deviations
+C for the event-to-event variation in the model parameters.
+C Note that these events are produced in the c.m. frame only.
+C
+C Anisotropic flow may also be simulated by introducing explicit
+C phi-dependence (azimuthal angle) in the particle distributions.
+C The assumed model is taken from Poskanzer and Voloshin, Phys. Rev.
+C C58, 1671 (1998), Eq.(1), where we use,
+C
+C E d^3N/dp^3 = (1/2*pi*pt)*[d^2N/dpt*dy]
+C * [1 + SUM(n=1,nflowterms){2*Vn*cos[n(phi-PSIr)]}]
+C
+C with up to 'nflowterms' (currently set to 6, see file
+C Parameter_values.inc) Fourier components allowed. Vn are
+C coefficients and PSIr is the reaction plane angle.
+C The algebraic signs of the Vn terms for n=odd are reversed
+C from their input values for particles with rapidity (y) < 0
+C as suggested in Poskanzer and Voloshin.
+C The flow parameters can depend on pt and rapidity (y) according
+C to the model suggested by Art Poskanzer (Feb. 2000) and as
+C defined in the Function Vn_pt_y.
+C
+C The user may select either to have the same multiplicity per
+C particle type for each event or to let the multiplicity vary
+C randomly according to a Poisson distribution. In addition, an
+C overall multiplicative scale factor can be applied to each
+C particle ID's multiplicity (same factor applied to each PID).
+C This scaling can vary randomly according to a Gaussian from
+C event-to-event. This is to simulate trigger acceptance
+C fluctuations. Similarly the
+C parameters of the one-particle distribution models may either
+C be fixed to the same value for each event or allowed to randomly
+C vary about a specified mean with a specified standard deviation
+C and assuming a Gaussian distribution.
+C
+C With respect to the reaction plane and anisotropic flow simulation,
+C the user may select among four options:
+C (1) ignore reaction plane and anisotropic flow effects; all
+C distributions will be azimuthally invariant, on average.
+C (2) assume a fixed reaction plane angle, PSIr, for all events
+C in the run.
+C (3) assume a Gaussian distribution with specified mean and
+C standard deviation for the reaction plane angles for the
+C events in the run. PSIr is randomly determined for each
+C event.
+C (4) assume uniformly distributed, random values for the reaction
+C plane angles from 0 to 360 deg., for each event in the run.
+C
+C The user may also select the anisotropic flow parameters, Vn,
+C to either be fixed for each event, or to randomly vary from event
+C to event according to a Gaussian distribution where the user must
+C specify the mean and std. dev. For both cases the input file must
+C list the 'nflowterms' (e.g. 6) values of the mean and Std.dev. for
+C the Vn parameters for all particle ID types included in the run.
+C
+C The available list of particles has been increased to permit a
+C number of meson and baryon resonances. For those with broad widths
+C the code samples the mass distribution for the resonance and outputs
+C the resonance mass for each instance in a special kinematic file
+C (see file unit=9, filename = 'mult_gen.kin'). The resonance shapes
+C are approximately Breit-Wigner and are specific for each resonance
+C case. The additional particle/resonances include: rho(+,-,0),
+C omega(0), eta', phi, J/Psi, Delta(-,0,+,++) and K*(+,-,0). Masses
+C are sampled for the rho, omega, phi, Deltas and D*s.
+C Refer to SUBR: Particle_prop and Particle_mass for the explicit
+C parameters, resonance shape models, and sampling ranges.
+C
+C The input is from a file, named 'mult_gen.in'. The output is
+C loaded into a file named 'mult_gen.out' which includes run
+C header information, event header information and the EVENT: and
+C TRACK: formats as in the new STAR TEXT Format for event generator
+C input to GSTAR. A log file, 'mult_gen.log' is also written which
+C may contain error messages. Normally this file should be empty
+C after a successful run. These filenames can easily be changed
+C to more suitable names by the script that runs the program or
+C by hand.
+C
+C The method for generating random multiplicities and model parameter
+C values involves the following steps:
+C (1) The Poisson or Gaussian distributions are computed and
+C loaded into function f().
+C (2) The distribution f(x') is integrated from xmin to x
+C and saved from x = xmin to x = xmax. The range and mesh
+C spaces are specified by the user.
+C (3) The integral of f is normalized to unity where
+C integral[f(x')](at x = xmin) = 0.0
+C integral[f(x')](at x = xmax) = 1.0
+C (4) A random number generator is called which delivers values
+C between 0.0 and 1.0.
+C (5) We consider the coordinate x (from xmin to xmax) to be
+C dependent on the integral[f]. Using the random number
+C for the selected value of integral[f] the value of x
+C is obtained by interpolation.
+C
+C An interpolation subroutine from Rubin Landau, Oregon State Univ.,
+C is used to do this interpolation; it involves uneven mesh point
+C spacing.
+C
+C The method for generating the particle momenta uses the
+C standard random elimination method and involves the following
+C steps:
+C
+C For model_type = 1,2,3,4 which are functions of pt,y (see following):
+C (1) The y range is computed using the pseudorapidity (eta)
+C range and includes ample cushioning around the sides
+C along the eta acceptance edges.
+C (2) The transverse momentum (pt) and rapidity (y) are
+C randomly chosen within the specified ranges.
+C (3) The pseudorapidity is computed for this (pt,y) value
+C (and the mass for each pid) and checked against the
+C pseudorapidity acceptance range.
+C (4) If the pseudorapidity is within range then the one-particle
+C model distribution is calculated at this point and its ratio
+C to the maximum value throughout (pt,eta) acceptance region
+C is calculated.
+C (5) Another random number is called and if less than the ratio
+C from step#4 the particle momentum is used; if not, then
+C another trial value of (pt,y) is obtained.
+C (6) This continues until the required multiplicity for the
+C specific event and particle type has been satisfied.
+C (7) This process is repeated for the requested number of particle
+C types and events.
+C
+C For model_type = 5,6 (see following) which are input bin-by-bin
+C in pt,eta:
+C (1) The transverse momentum (pt) and pseudorapidity (eta) are
+C randomly chosen within the specified ranges.
+C (2) The one-particle model distribution is calculated at this
+C point and its ratio to the maximum value throughout the
+C (pt,eta) region is calculated.
+C (3) Another random number is called and if less than the ratio
+C from step(2) the particle momentum is used; if not then
+C another trial value of (pt,eta) is obtained.
+C (4) This continues until the required multiplicity for the
+C specific event and particle type has been satisfied.
+C (5) This process is repeated for the requested number of particle
+C types and events.
+C
+C Problematic parameter values are tested, bad input values are checked
+C and in some cases may be changed so that the program will not crash.
+C In some cases the code execution is stopped.
+C Some distributions and/or unusual model parameter values may cause the
+C code to hang up due to the poor performance of the "elimination"
+C method for very strongly peaked distributions. These are tested for
+C certain problematic values and if necessary these events are aborted.
+C A message, "*** Event No. 2903 ABORTED:" for example is printed
+C in the 'mult_gen.out' file. Temperatures .le. 0.01 GeV and rapidity
+C width parameters .le. 0.01 will cause the event to abort.
+C
+C The input is described below in the 'read' statements and also in
+C the sample input file. Some additional comments are as follows:
+C
+C (1) n_events - Selected number of events in run. Can be anything
+C .ge. 1.
+C (2) n_pid_type - Number of particle ID types to include in the
+C particle list. e.g. pi(+) and pi(-) are counted
+C separately. The limit is set by parameter npid
+C in the accompanying include file 'Parameter_values.inc'
+C and is presently set at 20.
+C (3) model_type - equals 1,2,3,4,5 or 6 so far. See comments in
+C Function dNdpty to see what is calculated.
+C The models included are:
+C = 1, Factorized mt exponential, Gaussian rapidity model
+C = 2, Pratt non-expanding, spherical thermal source model
+C = 3, Bertsch non-expanding spherical thermal source model
+C = 4, Pratt spherically expanding, thermally equilibrated
+C source model.
+C = 5, Factorized pt and eta distributions input bin-by-bin.
+C = 6, Fully 2D pt,eta distributions input bin-by-bin.
+C NOTE: model_type = 1-4 are functions of (pt,y)
+C model_type = 5,6 are functions of (pt,eta)
+C (4) reac_plane_cntrl - Can be either 1,2,3 or 4 where:
+C = 1 to ignore reaction plane and anisotropic flow,
+C all distributions will be azimuthally symm.
+C = 2 to use a fixed reaction plane angle for all
+C events in the run.
+C = 3 to assume a randomly varying reaction plane
+C angle for each event as determined by a
+C Gaussian distribution.
+C = 4 to assume a randomly varying reaction plane
+C for each event in the run as determined by
+C a uniform distribution from 0 to 360 deg.
+C (5) PSIr_mean, PSIr_stdev - Reaction plane angle mean and Gaussian
+C std.dev. (both are in degrees) for cases
+C with reac_plane_cntrl = 2 (use mean value)
+C and 3. Note: these are read in regardless
+C of the value of reac_plane_cntrl.
+C (6) MultFac_mean, MultFac_stdev - Overall multiplicity scaling factor
+C for all PID types; mean and std.dev.;
+C for trigger fluctuations event-to-evt.
+C (7) pt_cut_min,pt_cut_max - Range of transverse momentum in GeV/c.
+C (8) eta_cut_min,eta_cut_max - Pseudorapidity range
+C (9) phi_cut_min,phi_cut_max - Azimuthal angular range in degrees.
+C (10) n_stdev_mult - Number of standard deviations about the mean value
+C of multiplicity to include in the random event-to-
+C event selection process. The maximum number of
+C steps that can be covered is determined by
+C parameter n_mult_max_steps in the accompanying
+C include file 'Parameter_values.inc' which is
+C presently set at 1000, but the true upper limit for
+C this is n_mult_max_steps - 1 = 999.
+C (11) n_stdev_temp - Same, except for the "Temperature" parameter.
+C (12) n_stdev_sigma- Same, except for the rapidity width parameter.
+C (13) n_stdev_expvel - Same, except for the expansion velocity parameter.
+C (14) n_stdev_PSIr - Same, except for the reaction plane angle
+C (15) n_stdev_Vn - Same, except for the anisotropy coefficients, Vn.
+C (16) n_stdev_MultFac - Same, except for the multiplicity scaling factor.
+C (17) n_integ_pts - Number of mesh points to use in the random model
+C parameter selection process. The upper limit is
+C set by parameter nmax_integ in the accompanying
+C include file 'Parameter_values.inc' which is presently
+C set at 100, but the true upper limit for n_integ_pts
+C is nmax_integ - 1 = 99.
+C (18) n_scan_pts - Number of mesh points to use to scan the (pt,y)
+C dependence of the model distributions looking for
+C the maximum value. The 2-D grid has
+C n_scan_pts * n_scan_pts points; no limit to size of
+C n_scan_pts.
+C (19) irand - Starting random number seed.
+C
+C**************************************************************************
+C FOR MODEL_TYPE = 1,2,3 or 4:
+C Input the following 7 lines for each particle type; repeat these
+C set of lines n_pid_type times:
+C
+C (a) gpid - Geant Particle ID code number
+C (b) mult_mean,mult_variance_control - Mean multiplicity and
+C variance control where:
+C mult_variance_control = 0 for no variance in multiplicity
+C mult_variance_control = 1 to allow Poisson distribution for
+C particle multiplicities for all events.
+C Note that a hard limit exists for the maximum possible
+C multiplicity for a given particle type per event. This is
+C determined by parameter factorial_max in accompanying include
+C file 'common_facfac.inc' and is presently set at 10000.
+C (c) Temp_mean, Temp_stdev - Temperature parameter mean (in GeV)
+C and standard deviation (Gaussian distribution assumed).
+C (d) sigma_mean, sigma_stdev - Rapidity distribution width (sigma)
+C parameter mean and standard deviation (Gaussian distribution
+C assumed).
+C (e) expvel_mean, expvel_stdev - S. Pratt expansion velocity
+C (in units of c) mean and standard deviation (Gaussian
+C distribution assumed).
+C (f) Vn_mean(k); k=1,4 - Anisotropic flow parameters, mean values
+C for Fourier component n=1.
+C (g) Vn_stdev(k); k=1,4 - Anisotropic flow parameters, std.dev.
+C values for Fourier component n=1.
+C
+C Repeat the last two lines of input for remaining Fourier
+C components n=2,3...6. Include all 6 sets of parameters
+C even if these are not used by the model for Vn(pt,y) (set
+C unused parameter means and std.dev. to 0.0). List 4 values
+C on every line, even though for n=even the 4th quantity is
+C not used.
+C
+C**************************************************************************
+C FOR MODEL_TYPE = 5 input the following set of lines for each particle
+C type; repeat these n_pid_type times.
+C
+C (a) gpid - Geant Particle ID code number
+C (b) mult_mean,mult_variance_control - Mean multiplicity and
+C variance control where:
+C mult_variance_control = 0 for no variance in multiplicity
+C mult_variance_control = 1 to allow Poisson distribution for
+C particle multiplicities for all events.
+C (c) pt_start, eta_start - minimum starting values for pt, eta
+C input for the bin-by-bin distributions.
+C (d) n_pt_bins, n_eta_bins - # input pt and eta bins.
+C (e) delta_pt, pt_bin - pt bin size and function value, repeat for
+C each pt bin.
+C (f) delta_eta, eta_bin - eta bin size and function value, repeat
+C for each eta bin.
+C (g) Vn_mean(k); k=1,4 - Anisotropic flow parameters, mean values
+C for Fourier component n=1.
+C (h) Vn_stdev(k); k=1,4 - Anisotropic flow parameters, std.dev.
+C values for Fourier component n=1.
+C
+C Repeat the last two lines of input for remaining Fourier
+C components n=2,3...6. Include all 6 sets of parameters
+C even if these are not used by the model for Vn(pt,y) (set
+C unused parameter means and std.dev. to 0.0). List 4 values
+C on every line, even though for n=even the 4th quantity is
+C not used.
+C
+C NOTE: The pt, eta ranges must fully include the requested ranges
+C in input #4 and 5 above; else the code execution will stop.
+C
+C Also, variable bin sizes are permitted for the input distributions.
+C
+C Also, this input distribution is used for all events in the run;
+C no fluctuations in this "parent" distribution are allowed from
+C event-to-event.
+C
+C**************************************************************************
+C FOR MODEL_TYPE = 6 input the following set of lines for each particle
+C type; repeat these n_pid_type times.
+C
+C (a) gpid - Geant Particle ID code number
+C (b) mult_mean,mult_variance_control - Mean multiplicity and
+C variance control where:
+C mult_variance_control = 0 for no variance in multiplicity
+C mult_variance_control = 1 to allow Poisson distribution for
+C particle multiplicities for all events.
+C (c) pt_start, eta_start - minimum starting values for pt, eta
+C input for the bin-by-bin distributions.
+C (d) n_pt_bins, n_eta_bins - # input pt and eta bins.
+C (e) delta_pt - pt bin size, repeat for each pt bin.
+C (f) delta_eta - eta bin size, repeat for each eta bin.
+C (g) i,j,pt_eta_bin(i,j) - read pt (index = i) and eta (index = j)
+C bin numbers and bin value for full 2D space.
+C (h) Vn_mean(k); k=1,4 - Anisotropic flow parameters, mean values
+C for Fourier component n=1.
+C (i) Vn_stdev(k); k=1,4 - Anisotropic flow parameters, std.dev.
+C values for Fourier component n=1.
+C
+C Repeat the last two lines of input for remaining Fourier
+C components n=2,3...6. Include all 6 sets of parameters
+C even if these are not used by the model for Vn(pt,y) (set
+C unused parameter means and std.dev. to 0.0). List 4 values
+C on every line, even though for n=even the 4th quantity is
+C not used.
+C
+C NOTE: The pt, eta ranges must fully include the requested ranges
+C in input #4 and 5 above; else the code execution will stop.
+C
+C Also, variable bin sizes are permitted for the input distributions.
+C
+C Also, this input distribution is used for all events in the run;
+C no fluctuations in this "parent" distribution are allowed from
+C event-to-event.
+C
+CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
+
+ Include 'Parameter_values.inc'
+ Include 'common_facfac.inc'
+ include 'common_dist_bin.inc'
+ Common/track/ pout(npid,4,factorial_max)
+ real*4 pout
+ Common/Geant/geant_mass(200),geant_charge(200),
+ 1 geant_lifetime(200),geant_width(200)
+ real*4 geant_mass,geant_charge,geant_lifetime,geant_width
+ Common/Mass/ Mass_integ_save(npid,nmax_integ),
+ 1 Mass_xfunc_save(npid,nmax_integ)
+ real*4 Mass_integ_save,Mass_xfunc_save
+
+ integer n_events, n_pid_type, model_type, n_integ_pts, irand
+ integer gpid(npid),mult_mean(npid),mult_variance_control(npid)
+ integer i,j,k,pid,mult_min,mult_max,n_mult_steps(npid),ievent
+ integer mult_event(npid),n_scan_pts,total_mult,n_vertices
+ integer event_abort,status_abort, status
+ integer iptbin, ietabin
+ integer track_counter,start_v,stop_v
+ integer reac_plane_cntrl
+ integer jodd,total_mult_inc(npid)
+
+ real*4 pt_cut_min,pt_cut_max,eta_cut_min,eta_cut_max
+ real*4 y_cut_min,y_cut_max
+ real*4 phi_cut_min,phi_cut_max,n_stdev_mult,n_stdev_temp
+ real*4 n_stdev_sigma,n_stdev_expvel,Temp_mean(npid)
+ real*4 Temp_stdev(npid),pi,rad,mult_mean_real,mult_stdev
+ real*4 mult_min_real,mult_max_real,ran
+ real*4 Temp_abort, sigma_abort, bin_value
+
+ real*4 mult_integ(n_mult_max_steps),mult_xfunc(n_mult_max_steps)
+ real*4 mult_integ_save(npid,n_mult_max_steps)
+ real*4 mult_xfunc_save(npid,n_mult_max_steps)
+ real*4 mult_event_real
+
+ real*4 Temp_min,Temp_max,integ(nmax_integ),xfunc(nmax_integ)
+ real*4 Temp_integ_save(npid,nmax_integ)
+ real*4 Temp_xfunc_save(npid,nmax_integ)
+ real*4 Temp_event(npid)
+
+ real*4 sigma_stdev(npid),sigma_mean(npid),sigma_min,sigma_max
+ real*4 sigma_integ_save(npid,nmax_integ)
+ real*4 sigma_xfunc_save(npid,nmax_integ)
+ real*4 sigma_event(npid)
+
+ real*4 expvel_stdev(npid), expvel_mean(npid)
+ real*4 expvel_min, expvel_max
+ real*4 expvel_integ_save(npid,nmax_integ)
+ real*4 expvel_xfunc_save(npid,nmax_integ)
+ real*4 expvel_event(npid)
+
+ real*4 masstemp,pxtemp,pytemp,pztemp,pttemp,etatemp,ytemp,phitemp
+
+CCC Variables associated with anisotropic flow:
+
+ real*4 PSIr_mean, PSIr_stdev, n_stdev_PSIr, n_stdev_Vn
+ real*4 PSIr_min, PSIr_max, PSIr_event
+ real*4 PSIr_integ_save(nmax_integ)
+ real*4 PSIr_xfunc_save(nmax_integ)
+ real*4 Vn_mean(nflowterms,4,npid), Vn_stdev(nflowterms,4,npid)
+ real*4 Vn_min, Vn_max
+ real*4 Vn1_integ_save(nflowterms,npid,nmax_integ)
+ real*4 Vn1_xfunc_save(nflowterms,npid,nmax_integ)
+ real*4 Vn2_integ_save(nflowterms,npid,nmax_integ)
+ real*4 Vn2_xfunc_save(nflowterms,npid,nmax_integ)
+ real*4 Vn3_integ_save(nflowterms,npid,nmax_integ)
+ real*4 Vn3_xfunc_save(nflowterms,npid,nmax_integ)
+ real*4 Vn4_integ_save(nflowterms,npid,nmax_integ)
+ real*4 Vn4_xfunc_save(nflowterms,npid,nmax_integ)
+ real*4 Vn_event(nflowterms,4,npid), Vn_event_tmp(nflowterms,4)
+ real*4 Vn_sum(nflowterms,npid),cosinefac(nflowterms)
+
+CCC Variables associated with trigger fluctuations:
+ real*4 MultFac_mean, MultFac_stdev, n_stdev_MultFac
+ real*4 MultFac_min, MultFac_max, MultFac_event
+ real*4 MultFac_integ_save(nmax_integ)
+ real*4 MultFac_xfunc_save(nmax_integ)
+CCC Open I/O Files:
+
+ open(unit=4,type='old',access='sequential',name='mult_gen.in')
+C--> Commented for AliRoot direct COMMON block access
+C open(unit=7,type='new',access='sequential',name='mult_gen.out')
+C open(unit=8,type='new',access='sequential',name='mult_gen.log')
+
+CCC NOTE:
+CCC Lines throughout the code which are commented out with "C-OUT"
+CCC can be activated to output special files (unit=9,10) with just the
+CCC mass,pt,eta,y,phi values listed for all the tracks for all events,
+CCC or the cos[n*(phi - PSI-reaction-plane)] for n=1,2...6.
+CCC These files can be directly loaded into PAW ntuples for analysis.
+
+C-OUT open(unit=9,type='new',access='sequential',name='mult_gen.kin')
+C-OUT open(unit=10,type='new',access='sequential',name='mult_gen.cos')
+C--> Commented for AliRoot direct COMMON block access
+C open(unit=9,type='new',access='sequential',name='mult_gen.kin')
+C open(unit=10,type='new',access='sequential',name='mult_gen.cos')
+
+CCC File 'mult_gen.in' is the input file for the run.
+CCC File 'mult_gen.out' is the output file in STAR New TEXT Format.
+CCC File 'mult_gen.log' is a log file for the run.
+CCC File 'mult_gen.kin', if activated, lists track kinematics for PAW ntuples.
+CCC File 'mult_gen.cos', if activated, lists the cos(n*phi) for PAW ntuples.
+
+CCC Initialize Arrays to Zero:
+
+ do i = 1,npid
+ gpid(i) = 0
+ mult_mean(i) = 0
+ mult_variance_control(i) = 0
+ n_mult_steps(i) = 0
+ Temp_mean(i) = 0.0
+ Temp_stdev(i) = 0.0
+ sigma_mean(i) = 0.0
+ sigma_stdev(i) = 0.0
+ expvel_mean(i) = 0.0
+ expvel_stdev(i) = 0.0
+ mult_event(i) = 0
+ Temp_event(i) = 0.0
+ sigma_event(i) = 0.0
+ expvel_event(i) = 0.0
+ total_mult_inc(i) = 0
+
+ do j = 1,n_mult_max_steps
+ mult_integ_save(i,j) = 0.0
+ mult_xfunc_save(i,j) = 0.0
+ end do
+
+ do j = 1,nmax_integ
+ Temp_integ_save(i,j) = 0.0
+ Temp_xfunc_save(i,j) = 0.0
+ sigma_integ_save(i,j) = 0.0
+ sigma_xfunc_save(i,j) = 0.0
+ expvel_integ_save(i,j) = 0.0
+ expvel_xfunc_save(i,j) = 0.0
+ Mass_integ_save(i,j) = 0.0
+ Mass_xfunc_save(i,j) = 0.0
+ end do
+ end do
+
+ do j = 1,nflowterms
+ cosinefac(j) = 0.0
+ do k = 1,4
+ Vn_event_tmp(j,k) = 0.0
+ end do
+ do i = 1,npid
+ Vn_sum(j,i) = 0.0
+ do k = 1,4
+ Vn_mean(j,k,i) = 0.0
+ Vn_stdev(j,k,i) = 0.0
+ Vn_event(j,k,i) = 0.0
+ end do
+ do k = 1,nmax_integ
+ Vn1_integ_save(j,i,k) = 0.0
+ Vn1_xfunc_save(j,i,k) = 0.0
+ Vn2_integ_save(j,i,k) = 0.0
+ Vn2_xfunc_save(j,i,k) = 0.0
+ Vn3_integ_save(j,i,k) = 0.0
+ Vn3_xfunc_save(j,i,k) = 0.0
+ Vn4_integ_save(j,i,k) = 0.0
+ Vn4_xfunc_save(j,i,k) = 0.0
+ end do
+ end do
+ end do
+
+ do i = 1,nmax_integ
+ PSIr_integ_save(i) = 0.0
+ PSIr_xfunc_save(i) = 0.0
+ MultFac_integ_save(i) = 0.0
+ MultFac_xfunc_save(i) = 0.0
+ end do
+
+ do i = 1,n_mult_max_steps
+ mult_integ(i) = 0.0
+ mult_xfunc(i) = 0.0
+ end do
+
+ do i = 1,nmax_integ
+ integ(i) = 0.0
+ xfunc(i) = 0.0
+ end do
+
+ do i = 1,factorial_max
+ FACLOG(i) = 0.0
+ end do
+
+ do i = 1,60
+ geant_mass(i) = 0.0
+ end do
+
+ do i = 1,npid
+ pt_start(i) = 0.0
+ pt_stop(i) = 0.0
+ eta_start(i) = 0.0
+ eta_stop(i) = 0.0
+ n_pt_bins(i) = 0
+ n_eta_bins(i) = 0
+ do j = 1,n_bins_max
+ delta_pt(i,j) = 0.0
+ delta_eta(i,j) = 0.0
+ pt_bin(i,j) = 0.0
+ eta_bin(i,j) = 0.0
+ do k = 1,n_bins_max
+ pt_eta_bin(i,j,k) = 0.0
+ end do
+ end do
+ end do
+
+CCC Read Input:
+
+ read(4,*) n_events ! No. of events to generate
+ read(4,*) n_pid_type ! No. of Geant PID types to include
+ read(4,*) model_type ! Distribution model type (see
+CCC ! Function dNdpty for explanation).
+ read(4,*) reac_plane_cntrl ! Reaction plane control option (1-4)
+ read(4,*) PSIr_mean,PSIr_stdev ! Reaction plane angle mean and std.
+CCC ! dev., both are in degrees.
+ read(4,*) MultFac_mean,MultFac_stdev ! Mult scaling factor mean,stdev.
+ read(4,*) pt_cut_min,pt_cut_max ! Min/Max pt range in GeV/c
+ read(4,*) eta_cut_min,eta_cut_max ! Min/Max pseudorapidity range
+ read(4,*) phi_cut_min,phi_cut_max ! Min/Max azimuthal angular range (deg)
+ read(4,*) n_stdev_mult ! No.(+/-) standard deviation range
+CCC ! for multiplicity
+ read(4,*) n_stdev_temp ! No.(+/-) st.dev. range for Temp.
+ read(4,*) n_stdev_sigma ! No.(+/-) st.dev. range for rapidity
+CCC ! width, sigma.
+ read(4,*) n_stdev_expvel ! No.(+/-) st.dev. range for expansion
+CCC ! velocity.
+ read(4,*) n_stdev_PSIr ! No.(+/-) st.dev. range for PSIr
+ read(4,*) n_stdev_Vn ! No.(+/-) st.dev. range for anisotropic
+CCC ! flow parameters Vn.
+ read(4,*) n_stdev_MultFac ! No.(+/-) st.dev. range for multipli-
+CCC ! city scaling factor for all PIDs.
+ read(4,*) n_integ_pts ! No. of integration mesh points to use
+CCC ! for random parameter fluctuations.
+ read(4,*) n_scan_pts ! No. of pt and eta mesh points to use
+CCC ! in scan for maximum value of dN/dpt*dy
+ read(4,*) irand ! Random number seed; default=12345.
+
+CCC Check Validity and Consistency of Input Parameters so far:
+
+ if(n_events .le. 0) n_events = 1
+ if(n_pid_type .le. 0) n_pid_type = 1
+ if(pt_cut_min .gt. pt_cut_max) then
+C--> write(8,40) pt_cut_min,pt_cut_max
+ RETURN
+ end if
+ if(eta_cut_min .gt. eta_cut_max) then
+C--> write(8,41) eta_cut_min,eta_cut_max
+ RETURN
+ end if
+ if(phi_cut_min .gt. phi_cut_max) then
+C--> write(8,42) phi_cut_min,phi_cut_max
+ RETURN
+ end if
+40 Format(//10x,'pt_cut_min = ',F7.4,' gt max = ',F7.4,' -STOP')
+41 Format(//10x,'eta_cut_min = ',F7.4,' gt max = ',F7.4,' -STOP')
+42 Format(//10x,'phi_cut_min = ',F7.4,' gt max = ',F7.4,' -STOP')
+ if(n_stdev_mult .lt. 0.0) n_stdev_mult = 1.0
+ if(n_stdev_temp .lt. 0.0) n_stdev_temp = 1.0
+ if(n_stdev_sigma .lt. 0.0) n_stdev_sigma = 1.0
+ if(n_stdev_expvel .lt. 0.0) n_stdev_expvel = 1.0
+ if(n_stdev_PSIr .lt. 0.0) n_stdev_PSIr = 1.0
+ if(n_stdev_Vn .lt. 0.0) n_stdev_Vn = 1.0
+ if(n_stdev_MultFac .lt. 0.0) n_stdev_MultFac = 1.0
+ if(n_integ_pts .le. 0) n_integ_pts = 10
+ if(n_scan_pts .le. 0) n_scan_pts = 10
+
+ if(irand .le. 0) irand = 12345
+ if(n_pid_type .gt. npid) then
+C--> write(8,10) n_pid_type, npid
+10 Format(//10x,'No. requested PID types = ',I7,
+ 1 ', exceeds maximum of ',I7,'; reset')
+ n_pid_type = npid
+ end if
+ if(model_type .lt. 0 .or. model_type .gt. 6) then
+C--> write(8,11) model_type
+11 Format(/10x,'model_type = ',I5,' is not allowed; STOP')
+ RETURN
+ end if
+ if(n_integ_pts .gt. nmax_integ) then
+C--> write(8,12) n_integ_pts, nmax_integ
+12 Format(/10x,'No. integ. pts = ',I7,
+ 1 ', exceeds maximum of ',I7,'; reset')
+ n_integ_pts = nmax_integ
+ end if
+ if(reac_plane_cntrl .lt. 1 .OR. reac_plane_cntrl .gt. 4) then
+C--> write(8,*) 'reac_plane_cntrl = ',reac_plane_cntrl,'-STOP'
+ RETURN
+ end if
+
+CCC Force the reaction plane angle (mean value) to be within the
+CCC range 0 to 360 deg.
+ PSIr_mean = PSIr_mean - 360.0*float(int(PSIr_mean/360.0))
+ if(PSIr_mean .lt. 0.0) PSIr_mean = PSIr_mean + 360.0
+ if(PSIr_stdev .lt. 0.0) then
+C--> write(8,*) 'Reaction plane angle Std. dev. is < 0.0 - STOP'
+ RETURN
+ end if
+
+CCC Check the multiplicity scaling factor input parameters:
+ if(MultFac_mean.lt.0.0 .or. MultFac_stdev.lt.0.0) then
+C--> write(8,48) MultFac_mean, MultFac_stdev
+48 Format(//10x,'Multiplicity scaling factor mean or stdev = ',
+ 1 2F7.4,' - not valid - STOP')
+ RETURN
+ end if
+
+CCC FOR MODEL_TYPE = 1,2,3 or 4;
+CCC Repeat the following lines of input for each particle ID type:
+
+ IF(model_type.ge.1 .and. model_type.le.4) then
+
+ do pid = 1,n_pid_type
+
+ read(4,*) gpid(pid) ! Geant Particle ID code number
+
+ read(4,*) mult_mean(pid), mult_variance_control(pid)
+CCC Mean multiplicity and variance control where:
+CCC mult_variance_control = 0 for no variance in multiplicity
+CCC mult_variance_control = 1 to allow Poisson distribution for
+CCC particle multiplicities for all events.
+
+ read(4,*) Temp_mean(pid), Temp_stdev(pid)
+CCC Temperature parameter mean (in GeV) and standard deviation (Gaussian
+CCC distribution assumed).
+
+ read(4,*) sigma_mean(pid), sigma_stdev(pid)
+CCC Rapidity distribution width (sigma) parameter mean and standard
+CCC deviation (Gaussian distribution assumed).
+
+ read(4,*) expvel_mean(pid), expvel_stdev(pid)
+CCC S. Pratt expansion velocity (in units of c) mean and standard
+CCC deviation (Gaussian distribution assumed).
+
+ do i = 1,nflowterms
+ read(4,*) (Vn_mean(i,k,pid) ,k=1,4)
+ read(4,*) (Vn_stdev(i,k,pid),k=1,4)
+ end do
+CCC Anisotropic flow Fourier component coefficients specifying the
+CCC azimuthal (phi angle) dependence as defined in Phys. Rev. C58,
+CCC 1671 (1998), Eq.(1); 'nflowterms' number of parameters are
+CCC allowed (see file 'Parameter_values.inc' where this is currently
+CCC set to 6); these are the mean and Gaussian std. dev. to be used
+CCC if random, Gaussian variation in the Vn values from event-to-event
+CCC are selected (via reac_plane_cntrl).
+CCC The flow parameters are pt and y dependent; see Function Vn_pt_y
+CCC for the full definitions.
+
+CCC Check Validity and Consistency of Input Parameters
+
+ if(Temp_mean(pid).le.0.0 .or. Temp_stdev(pid).lt.0.0) then
+C--> write(8,45) pid,Temp_mean(pid),Temp_stdev(pid)
+45 Format(//10x,'For pid # ',I7,', Temp_mean or Temp_stdev= ',
+ 1 2F7.4,' - not valid -STOP')
+ RETURN
+ end if
+ if(sigma_mean(pid).le.0.0 .or. sigma_stdev(pid).lt.0.0) then
+C--> write(8,46) pid,sigma_mean(pid),sigma_stdev(pid)
+46 Format(//10x,'For pid # ',I7,', sigma_mean or sigma_stdev= ',
+ 1 2F7.4,' - not valid -STOP')
+ RETURN
+ end if
+ if(expvel_mean(pid).lt.0.0.or.expvel_stdev(pid).lt.0.0) then
+C--> write(8,47) pid,expvel_mean(pid),expvel_stdev(pid)
+47 Format(//10x,'For pid #',I7,',expvel_mean or expvel_stdev=',
+ 1 2F7.4,' - not valid -STOP')
+ RETURN
+ end if
+
+ do k = 1,4
+ do i = 1,nflowterms
+ if(Vn_stdev(i,k,pid) .lt. 0.0) then
+C--> write(8,*) 'For pid#,term#,k= ',pid,i,k
+C--> write(8,*) 'Vn_stdev = ',Vn_stdev(i,k,pid)
+C--> write(8,*) 'which is not valid, must be => 0, STOP'
+ RETURN
+ end if
+ end do
+ end do
+
+ end do ! End of loop over Particle ID types.
+
+ ELSE IF (model_type .eq. 5) then
+
+CCC Input for Factorized pt, eta bin-by-bin distribution:
+
+ do pid = 1,n_pid_type
+ read(4,*) gpid(pid)
+ read(4,*) mult_mean(pid), mult_variance_control(pid)
+ read(4,*) pt_start(pid), eta_start(pid)
+ read(4,*) n_pt_bins(pid), n_eta_bins(pid)
+ do i = 1,n_pt_bins(pid)
+ read(4,*) delta_pt(pid,i), pt_bin(pid,i)
+ end do
+ do i = 1,n_eta_bins(pid)
+ read(4,*) delta_eta(pid,i), eta_bin(pid,i)
+ end do
+
+ do i = 1,nflowterms
+ read(4,*) (Vn_mean(i,k,pid) ,k=1,4)
+ read(4,*) (Vn_stdev(i,k,pid),k=1,4)
+ end do
+CCC Anisotropic flow Fourier component coefficients specifying the
+CCC azimuthal (phi angle) dependence as defined in Phys. Rev. C58,
+CCC 1671 (1998), Eq.(1); 'nflowterms' number of parameters are
+CCC allowed (see file 'Parameter_values.inc' where this is currently
+CCC set to 6); these are the mean and Gaussian std. dev. to be used
+CCC if random, Gaussian variation in the Vn values from event-to-event
+CCC are selected (via reac_plane_cntrl).
+CCC The flow parameters are pt and y dependent; see Function Vn_pt_y
+CCC for the full definitions.
+
+CCC Evaluate grid and find maximum values of pt and eta for input bins:
+
+ pt_stop(pid) = pt_start(pid)
+ do i = 1,n_pt_bins(pid)
+ pt_stop(pid) = pt_stop(pid) + delta_pt(pid,i)
+ pt_bin_mesh(pid,i) = pt_stop(pid)
+ end do
+ eta_stop(pid) = eta_start(pid)
+ do i = 1,n_eta_bins(pid)
+ eta_stop(pid) = eta_stop(pid) + delta_eta(pid,i)
+ eta_bin_mesh(pid,i) = eta_stop(pid)
+ end do
+
+CCC Check ranges of pt and eta coverage:
+
+ if(pt_cut_min .lt. pt_start(pid)) then
+C--> write(8,50) pt_cut_min,pt_start(pid)
+50 Format(//10x,'pt_cut_min = ',F10.7,' .lt. pt_start =',
+ 1 F10.7,' - STOP')
+ RETURN
+ end if
+ if(pt_cut_max .gt. pt_stop(pid)) then
+C--> write(8,51) pt_cut_max,pt_stop(pid)
+51 Format(//10x,'pt_cut_max = ',F10.7,' .gt. pt_stop =',
+ 1 F10.7,' - STOP')
+ RETURN
+ end if
+ if(eta_cut_min .lt. eta_start(pid)) then
+C--> write(8,52) eta_cut_min,eta_start(pid)
+52 Format(//10x,'eta_cut_min = ',F10.7,'.lt.eta_start =',
+ 1 F10.7,' - STOP')
+ RETURN
+ end if
+ if(eta_cut_max .gt. eta_stop(pid)) then
+C--> write(8,53) eta_cut_max,eta_stop(pid)
+53 Format(//10x,'eta_cut_max = ',F10.7,'.gt.eta_stop =',
+ 1 F10.7,' - STOP')
+ RETURN
+ end if
+
+ do k = 1,4
+ do i = 1,nflowterms
+ if(Vn_stdev(i,k,pid) .lt. 0.0) then
+C--> write(8,*) 'For pid#,term#,k= ',pid,i,k
+C--> write(8,*) 'Vn_stdev = ',Vn_stdev(i,k,pid)
+C--> write(8,*) 'which is not valid, must be => 0, STOP'
+ RETURN
+ end if
+ end do
+ end do
+
+ end do ! End loop over particle ID types.
+
+ ELSE IF (model_type .eq. 6) then
+
+CCC Input for Full 2D (pt,eta) bin-by-bin distribution:
+
+ do pid = 1,n_pid_type
+ read(4,*) gpid(pid)
+ read(4,*) mult_mean(pid), mult_variance_control(pid)
+ read(4,*) pt_start(pid), eta_start(pid)
+ read(4,*) n_pt_bins(pid), n_eta_bins(pid)
+ do i = 1,n_pt_bins(pid)
+ read(4,*) delta_pt(pid,i)
+ end do
+ do i = 1,n_eta_bins(pid)
+ read(4,*) delta_eta(pid,i)
+ end do
+
+CCC Evaluate grid and find maximum values of pt and eta for input bins:
+
+ pt_stop(pid) = pt_start(pid)
+ do i = 1,n_pt_bins(pid)
+ pt_stop(pid) = pt_stop(pid) + delta_pt(pid,i)
+ pt_bin_mesh(pid,i) = pt_stop(pid)
+ end do
+ eta_stop(pid) = eta_start(pid)
+ do i = 1,n_eta_bins(pid)
+ eta_stop(pid) = eta_stop(pid) + delta_eta(pid,i)
+ eta_bin_mesh(pid,i) = eta_stop(pid)
+ end do
+
+CCC Load 2D bin-by-bin distribution:
+
+ do i = 1,n_pt_bins(pid)*n_eta_bins(pid)
+ read(4,*) iptbin,ietabin,bin_value
+ pt_eta_bin(pid,iptbin,ietabin) = bin_value
+ end do
+
+ do i = 1,nflowterms
+ read(4,*) (Vn_mean(i,k,pid) ,k=1,4)
+ read(4,*) (Vn_stdev(i,k,pid),k=1,4)
+ end do
+CCC Anisotropic flow Fourier component coefficients specifying the
+CCC azimuthal (phi angle) dependence as defined in Phys. Rev. C58,
+CCC 1671 (1998), Eq.(1); 'nflowterms' number of parameters are
+CCC allowed (see file 'Parameter_values.inc' where this is currently
+CCC set to 6); these are the mean and Gaussian std. dev. to be used
+CCC if random, Gaussian variation in the Vn values from event-to-event
+CCC are selected (via reac_plane_cntrl).
+CCC The flow parameters are pt and y dependent; see Function Vn_pt_y
+CCC for the full definitions.
+
+CCC Check ranges of pt and eta coverage:
+
+ if(pt_cut_min .lt. pt_start(pid)) then
+C--> write(8,50) pt_cut_min,pt_start(pid)
+ RETURN
+ end if
+ if(pt_cut_max .gt. pt_stop(pid)) then
+C--> write(8,51) pt_cut_max,pt_stop(pid)
+ RETURN
+ end if
+ if(eta_cut_min .lt. eta_start(pid)) then
+C--> write(8,52) eta_cut_min,eta_start(pid)
+ RETURN
+ end if
+ if(eta_cut_max .gt. eta_stop(pid)) then
+C--> write(8,53) eta_cut_max,eta_stop(pid)
+ RETURN
+ end if
+ do k = 1,4
+ do i = 1,nflowterms
+ if(Vn_stdev(i,k,pid) .lt. 0.0) then
+C--> write(8,*) 'For pid#,term#,k= ',pid,i,k
+C--> write(8,*) 'Vn_stdev = ',Vn_stdev(i,k,pid)
+C--> write(8,*) 'which is not valid, must be => 0, STOP'
+ RETURN
+ end if
+ end do
+ end do
+
+ end do ! End loop over particle ID types.
+
+ END IF ! End of MODEL_TYPE Options input:
+
+CCC Check some more input parameters; Set constants, and load data tables:
+
+ do pid = 1,n_pid_type
+ if(gpid(pid).le.0 .or. gpid(pid).gt.200) then
+C--> write(8,43) pid,gpid(pid)
+43 Format(//10x,'For pid # ',I7,', gpid = ',I7,' - not valid -STOP')
+ RETURN
+ end if
+ if(mult_variance_control(pid) .lt. 0 .or.
+ 1 mult_variance_control(pid) .gt. 1)
+ 2 mult_variance_control(pid) = 0
+ if(mult_mean(pid) .le. 0) then
+C--> write(8,44) pid,mult_mean(pid)
+44 Format(//10x,'For pid # ',I7,', mult_mean= ',I7,
+ 1 ' - not valid -STOP')
+ RETURN
+ end if
+ end do ! End Particle ID input parameter check and verification loop.
+
+ pi = 3.141592654
+ rad = 180.0/pi
+ Temp_abort = 0.01
+ sigma_abort = 0.01
+
+CCC Load particle properties array; mass in GeV, etc.:
+
+ Call Particle_prop
+
+CCC Load log(n!) values into Common/FACFAC/
+
+ Call FACTORIAL
+
+CCC Set y (rapidity) range, to be used for model_type = 1-4
+ if(eta_cut_min .ge. 0.0) then
+ y_cut_min = 0.0
+ else
+ y_cut_min = eta_cut_min
+ end if
+ if(eta_cut_max .le. 0.0) then
+ y_cut_max = 0.0
+ else
+ y_cut_max = eta_cut_max
+ end if
+
+CCC Obtain integrals of 1D distribution functions of multiplicity
+CCC (Poisson, integer) and parameters (Gaussian, real*4) for the
+CCC particle model distributions, for each particle ID type.
+CCC These integrated quantities are then used with random number
+CCC selection to generate a distribution of multiplicities and
+CCC parameter values for each event in the run.
+
+CCC Obtain 1D integrals for Gaussian distributions for reaction plane
+CCC angles:
+
+ if(reac_plane_cntrl .eq. 3) then
+ if((n_stdev_PSIr*PSIr_stdev).gt. 0.0) then
+ PSIr_min = PSIr_mean - n_stdev_PSIr*PSIr_stdev
+ PSIr_max = PSIr_mean + n_stdev_PSIr*PSIr_stdev
+ Call Gaussian(PSIr_min,PSIr_max,PSIr_mean,PSIr_stdev,
+ 1 n_integ_pts,integ,xfunc,nmax_integ)
+ do i = 1,n_integ_pts + 1
+ PSIr_integ_save(i) = integ(i)
+ PSIr_xfunc_save(i) = xfunc(i)
+ end do
+ else
+ PSIr_event = PSIr_mean
+ end if
+ end if
+
+CCC Obtain 1D integral for Gaussian distribution for the trigger fluctuation
+CCC simulations via the overall multiplicity scaling factor.
+ if((n_stdev_MultFac*MultFac_stdev).gt.0.0) then
+ Call MinMax(MultFac_mean,MultFac_stdev,n_stdev_MultFac,
+ 1 MultFac_min,MultFac_max)
+ Call Gaussian(MultFac_min,MultFac_max,MultFac_mean,
+ 1 MultFac_stdev,n_integ_pts,integ,xfunc,nmax_integ)
+ do i = 1,n_integ_pts + 1
+ MultFac_integ_save(i) = integ(i)
+ MultFac_xfunc_save(i) = xfunc(i)
+ end do
+ else
+ MultFac_event = MultFac_mean
+ end if
+
+ do pid = 1,n_pid_type
+
+ Call Particle_mass(gpid(pid),pid,n_integ_pts)
+
+ if(mult_variance_control(pid) .ne. 0) then
+ mult_mean_real = float(mult_mean(pid))
+ mult_stdev = sqrt(mult_mean_real)
+ Call MinMax(mult_mean_real,mult_stdev,n_stdev_mult,
+ 1 mult_min_real,mult_max_real)
+ mult_min = int(mult_min_real)
+ mult_max = int(mult_max_real + 1)
+ n_mult_steps(pid) = mult_max - mult_min + 1
+ if((n_mult_steps(pid) + 1) .gt. n_mult_max_steps) then
+C--> write(8,20) n_mult_steps(pid),n_mult_max_steps
+20 Format(//10x,'No. steps in multiplicity integral = ',
+ 1 I7,' + 1, exceeds max no. of ',I7,'; reset')
+ mult_min = mult_mean(pid) - (n_mult_max_steps - 1)/2
+ mult_max = mult_mean(pid) + (n_mult_max_steps - 1)/2
+ n_mult_steps(pid) = mult_max - mult_min + 1
+ end if
+ if((mult_max + 1) .gt. factorial_max) then
+C--> write(8,30) mult_max, factorial_max
+30 Format(//10x,'In Poisson multiplicity distribution,',
+ 1 ' max = ',I7,', exceeds limit of ',I7,' - STOP')
+ RETURN
+ end if
+
+ Call Poisson(mult_min,mult_max,mult_mean(pid),
+ 1 n_mult_steps(pid),mult_integ,mult_xfunc,n_mult_max_steps)
+ do i = 1,n_mult_steps(pid) + 1
+ mult_integ_save(pid,i) = mult_integ(i)
+ mult_xfunc_save(pid,i) = mult_xfunc(i)
+ end do
+ else if (mult_variance_control(pid) .eq. 0) then
+ mult_event(pid) = mult_mean(pid)
+ end if
+
+ if(model_type .le. 4) then
+ if(Temp_stdev(pid) .ne. 0.0) then
+ Call MinMax(Temp_mean(pid),Temp_stdev(pid),n_stdev_temp,
+ 1 Temp_min, Temp_max)
+ Call Gaussian(Temp_min,Temp_max,Temp_mean(pid),
+ 1 Temp_stdev(pid),n_integ_pts,integ,xfunc,nmax_integ)
+ do i = 1,n_integ_pts + 1
+ Temp_integ_save(pid,i) = integ(i)
+ Temp_xfunc_save(pid,i) = xfunc(i)
+ end do
+ else if(Temp_stdev(pid) .eq. 0.0) then
+ Temp_event(pid) = Temp_mean(pid)
+ end if
+
+ if(sigma_stdev(pid) .ne. 0.0) then
+ Call MinMax(sigma_mean(pid),sigma_stdev(pid),n_stdev_sigma,
+ 1 sigma_min, sigma_max)
+ Call Gaussian(sigma_min,sigma_max,sigma_mean(pid),
+ 1 sigma_stdev(pid),n_integ_pts,integ,xfunc,nmax_integ)
+ do i = 1,n_integ_pts + 1
+ sigma_integ_save(pid,i) = integ(i)
+ sigma_xfunc_save(pid,i) = xfunc(i)
+ end do
+ else if(sigma_stdev(pid) .eq. 0.0) then
+ sigma_event(pid) = sigma_mean(pid)
+ end if
+
+ if(expvel_stdev(pid) .ne. 0.0) then
+ Call MinMax(expvel_mean(pid),expvel_stdev(pid),n_stdev_expvel,
+ 1 expvel_min, expvel_max)
+ Call Gaussian(expvel_min,expvel_max,expvel_mean(pid),
+ 1 expvel_stdev(pid),n_integ_pts,integ,xfunc,nmax_integ)
+ do i = 1,n_integ_pts + 1
+ expvel_integ_save(pid,i) = integ(i)
+ expvel_xfunc_save(pid,i) = xfunc(i)
+ end do
+ else if(expvel_stdev(pid) .eq. 0.0) then
+ expvel_event(pid) = expvel_mean(pid)
+ end if
+ end if ! End model_type .le. 4 options.
+
+ if(reac_plane_cntrl .gt. 1) then
+ do i = 1,nflowterms
+ do k = 1,4
+ if((n_stdev_Vn*Vn_stdev(i,k,pid)) .gt. 0.0) then
+ Vn_min = Vn_mean(i,k,pid)-n_stdev_Vn*Vn_stdev(i,k,pid)
+ Vn_max = Vn_mean(i,k,pid)+n_stdev_Vn*Vn_stdev(i,k,pid)
+ Call Gaussian(Vn_min,Vn_max,Vn_mean(i,k,pid),
+ 1 Vn_stdev(i,k,pid),n_integ_pts,integ,xfunc,nmax_integ)
+ if(k.eq.1) then
+ do j = 1,n_integ_pts + 1
+ Vn1_integ_save(i,pid,j) = integ(j)
+ Vn1_xfunc_save(i,pid,j) = xfunc(j)
+ end do
+ else if(k.eq.2) then
+ do j = 1,n_integ_pts + 1
+ Vn2_integ_save(i,pid,j) = integ(j)
+ Vn2_xfunc_save(i,pid,j) = xfunc(j)
+ end do
+ else if(k.eq.3) then
+ do j = 1,n_integ_pts + 1
+ Vn3_integ_save(i,pid,j) = integ(j)
+ Vn3_xfunc_save(i,pid,j) = xfunc(j)
+ end do
+ else if(k.eq.4) then
+ do j = 1,n_integ_pts + 1
+ Vn4_integ_save(i,pid,j) = integ(j)
+ Vn4_xfunc_save(i,pid,j) = xfunc(j)
+ end do
+ end if
+ else
+ Vn_event(i,k,pid) = Vn_mean(i,k,pid)
+ end if
+ end do
+ end do
+ end if
+
+ end do ! End of PID Loop:
+
+CCC Write Run Header Output:
+
+C--> write(7,200)
+C--> write(7,201) n_events,n_pid_type
+C--> if(model_type .eq. 1) write(7,202)
+C--> if(model_type .eq. 2) write(7,203)
+C--> if(model_type .eq. 3) write(7,204)
+C--> if(model_type .eq. 4) write(7,205)
+C--> if(model_type .eq. 5) write(7,2051)
+C--> if(model_type .eq. 6) write(7,2052)
+C--> write(7,2053) reac_plane_cntrl
+C--> write(7,2054) PSIr_mean, PSIr_stdev
+C--> write(7,2055) MultFac_mean,MultFac_stdev
+C--> write(7,206) pt_cut_min, pt_cut_max
+C--> write(7,207) eta_cut_min, eta_cut_max
+C--> write(7,2071) y_cut_min,y_cut_max
+C--> write(7,208) phi_cut_min, phi_cut_max
+C--> write(7,209) n_stdev_mult,n_stdev_temp,n_stdev_sigma,
+C--> 1 n_stdev_expvel
+C--> write(7,2091) n_stdev_PSIr, n_stdev_Vn
+C--> write(7,2092) n_stdev_MultFac
+C--> write(7,210) n_integ_pts
+C--> write(7,211) n_scan_pts
+C--> write(7,212) irand
+C--> write(7,213) (gpid(pid),pid = 1,n_pid_type)
+C--> write(7,214) (mult_mean(pid),pid = 1,n_pid_type)
+C--> write(7,215) (mult_variance_control(pid),pid = 1,n_pid_type)
+C--> if(model_type .le. 4) then
+C--> write(7,216) (Temp_mean(pid),pid = 1,n_pid_type)
+C--> write(7,217) (Temp_stdev(pid),pid = 1,n_pid_type)
+C--> write(7,218) (sigma_mean(pid),pid = 1,n_pid_type)
+C--> write(7,219) (sigma_stdev(pid),pid = 1,n_pid_type)
+C--> write(7,220) (expvel_mean(pid),pid = 1,n_pid_type)
+C--> write(7,221) (expvel_stdev(pid),pid = 1,n_pid_type)
+C--> else if(model_type .ge. 5) then
+C--> write(7,222) (n_pt_bins(pid),pid = 1,n_pid_type)
+C--> write(7,223) (n_eta_bins(pid),pid = 1,n_pid_type)
+C--> write(7,224) (pt_start(pid),pid = 1,n_pid_type)
+C--> write(7,225) (pt_stop(pid),pid = 1,n_pid_type)
+C--> write(7,226) (eta_start(pid),pid = 1,n_pid_type)
+C--> write(7,227) (eta_stop(pid),pid = 1,n_pid_type)
+C--> end if
+CCC Print out flow parameters:
+ do pid = 1,n_pid_type
+ do i = 1,nflowterms
+C--> write(7,2271) pid,i,(Vn_mean(i,k,pid),k=1,4)
+C--> write(7,2272) pid,i,(Vn_stdev(i,k,pid),k=1,4)
+ end do
+ end do
+
+200 Format('*** Multiplicity Generator Run Header ***')
+201 Format('* Number of events = ',I7,'; number of Particle ID',
+ 1 ' types = ',I5)
+202 Format('* Factorized mt exponential, Gaussian rapidity model')
+203 Format('* Pratt non-expanding, spherical thermal source model')
+204 Format('* Bertsch non-expanding spherical thermal source model')
+205 Format('* Pratt spherically expanding, thermally equilibrated ',
+ 1 'source model')
+2051 Format('* Factorized pt,eta bin-by-bin Distribution')
+2052 Format('* Full 2D pt,eta bin-by-bin Distribution')
+2053 Format('* Reaction plane control = ',I5)
+2054 Format('* Reaction plane angles - mean and std.dev. = ',2G12.5,
+ 1 ' (deg)')
+2055 Format('* Multiplicity Scaling Factor - mean and std.dev = ',
+ 1 2G12.5)
+206 Format('* Min, Max range in pt = ', 2G12.5)
+207 Format('* Min, Max range in pseudorapidity = ', 2G12.5)
+2071 Format('* Min, Max range in rapdity + cush = ', 2G12.5)
+208 Format('* Min, Max range in azimuthal angle = ', 2G12.5)
+209 Format('* No. std. dev. range used for mult and parameters = ',
+ 1 4F5.2)
+2091 Format('* No. std. dev. range for PSIr, Vn = ', 2G12.5)
+2092 Format('* No. std. dev. range for MultFac = ',G12.5)
+210 Format('* No. integration points for parameter variance = ',
+ 1 I6)
+211 Format('* No. of dN/dp(pt,y) scanning grid points to find ',
+ 1 'maximum = ', I6)
+212 Format('* Starting Random Number Seed = ',I10)
+213 Format('* Geant PID: ',10I7)
+214 Format('* Mean Multiplicity: ',10I7)
+215 Format('* Mult. Var. Control: ',10I7)
+216 Format('* Mean Temperature: ',10F7.4)
+217 Format('* Std. Dev. Temp: ',10F7.4)
+218 Format('* Mean Rap. sigma: ',10F7.4)
+219 Format('* Std. Dev. y-sigma: ',10F7.4)
+220 Format('* Mean expansion vel: ',10F7.4)
+221 Format('* Std. Dev. exp. vel: ',10F7.4)
+222 Format('* No. input pt bins: ',10I7)
+223 Format('* No. input eta bins: ',10I7)
+224 Format('* Input pt start: ',10F7.4)
+225 Format('* Input pt stop: ',10F7.4)
+226 Format('* Input eta start: ',10F7.4)
+227 Format('* Input eta stop: ',10F7.4)
+2271 Format('* Anisotropy Vn mean (pid=',I2,',n=',I1,'):',4G12.5)
+2272 Format('* Anisotropy Vn std. (pid=',I2,',n=',I1,'):',4G12.5)
+
+CCC END Run Header Output
+
+C**************************************************************
+C **
+C START EVENT LOOP **
+C **
+C**************************************************************
+
+ DO ievent = 1,n_events
+
+CCC Compute the Reaction plane angle for this event:
+ if(reac_plane_cntrl .eq. 1) then
+ PSIr_event = 0.0
+ else if(reac_plane_cntrl .eq. 2) then
+ PSIr_event = PSIr_mean
+ else if(reac_plane_cntrl .eq. 3) then
+ if((n_stdev_PSIr*PSIr_stdev) .gt. 0.0) then
+ do i = 1,n_integ_pts + 1
+ integ(i) = PSIr_integ_save(i)
+ xfunc(i) = PSIr_xfunc_save(i)
+ end do
+ Call LAGRNG(ran(irand),integ,PSIr_event,xfunc,
+ 1 n_integ_pts+1,1,5,n_integ_pts+1,1)
+CCC Ensure that the randomly selected reaction plane angle
+CCC is within the 0 to 360 deg range.
+ PSIr_event=PSIr_event-360.0*float(int(PSIr_event/360.0))
+ if(PSIr_event .lt. 0.0) PSIr_event = PSIr_event + 360.0
+ end if
+CCC NOTE: If PSIr_stdev=0.0, PSIr_event already calculated (see preceding)
+ else if(reac_plane_cntrl .eq. 4) then
+ PSIr_event = 360.0*ran(irand)
+ else
+ PSIr_event = 0.0
+ end if
+
+CCC Compute the multiplicity scaling factor for simulating trigger
+CCC fluctuations for this event:
+ if((n_stdev_MultFac*MultFac_stdev).gt.0.0) then
+ do i = 1,n_integ_pts + 1
+ integ(i) = MultFac_integ_save(i)
+ xfunc(i) = MultFac_xfunc_save(i)
+ end do
+ Call LAGRNG(ran(irand),integ,MultFac_event,xfunc,
+ 1 n_integ_pts+1,1,5,n_integ_pts+1,1)
+ end if
+
+ do pid = 1,n_pid_type
+ if(mult_variance_control(pid) .ne. 0) then
+ do i = 1,n_mult_steps(pid) + 1
+ mult_integ(i) = mult_integ_save(pid,i)
+ mult_xfunc(i) = mult_xfunc_save(pid,i)
+ end do
+ Call LAGRNG(ran(irand),mult_integ,mult_event_real,
+ 1 mult_xfunc,n_mult_steps(pid)+1,1,5,
+ 2 n_mult_steps(pid)+1,1)
+ mult_event(pid) = mult_event_real
+ else if(mult_variance_control(pid) .eq. 0) then
+ mult_event(pid) = mult_mean(pid)
+ end if
+ mult_event(pid) = int(MultFac_event*float(mult_event(pid))
+ 1 + 0.5)
+CCC Check each multiplicity wrt upper array size limit:
+ if(mult_event(pid).gt.factorial_max) then
+C--> write(8,*) 'For event#',ievent,'and pid#',pid,
+C--> 1 'multiplicity=',mult_event(pid),
+C--> 2 'which is > array size (factorial_max=',
+C--> 3 factorial_max,')-STOP'
+ RETURN
+ end if
+
+ if(model_type .le. 4) then
+
+ if(Temp_stdev(pid) .ne. 0.0) then
+ do i = 1,n_integ_pts + 1
+ integ(i) = Temp_integ_save(pid,i)
+ xfunc(i) = Temp_xfunc_save(pid,i)
+ end do
+ Call LAGRNG(ran(irand),integ,Temp_event(pid),xfunc,
+ 1 n_integ_pts+1,1,5,n_integ_pts+1,1)
+ end if
+
+ if(sigma_stdev(pid) .ne. 0.0) then
+ do i = 1,n_integ_pts + 1
+ integ(i) = sigma_integ_save(pid,i)
+ xfunc(i) = sigma_xfunc_save(pid,i)
+ end do
+ Call LAGRNG(ran(irand),integ,sigma_event(pid),xfunc,
+ 1 n_integ_pts+1,1,5,n_integ_pts+1,1)
+ end if
+
+ if(expvel_stdev(pid) .ne. 0.0) then
+ do i = 1,n_integ_pts + 1
+ integ(i) = expvel_integ_save(pid,i)
+ xfunc(i) = expvel_xfunc_save(pid,i)
+ end do
+ Call LAGRNG(ran(irand),integ,expvel_event(pid),xfunc,
+ 1 n_integ_pts+1,1,5,n_integ_pts+1,1)
+ end if
+ end if
+
+ if(reac_plane_cntrl .gt. 1) then
+ do i = 1,nflowterms
+ do k = 1,4
+ if((n_stdev_Vn*Vn_stdev(i,k,pid)) .gt. 0.0) then
+ if(k.eq.1) then
+ do j = 1,n_integ_pts + 1
+ integ(j) = Vn1_integ_save(i,pid,j)
+ xfunc(j) = Vn1_xfunc_save(i,pid,j)
+ end do
+ else if (k.eq.2) then
+ do j = 1,n_integ_pts + 1
+ integ(j) = Vn2_integ_save(i,pid,j)
+ xfunc(j) = Vn2_xfunc_save(i,pid,j)
+ end do
+ else if (k.eq.3) then
+ do j = 1,n_integ_pts + 1
+ integ(j) = Vn3_integ_save(i,pid,j)
+ xfunc(j) = Vn3_xfunc_save(i,pid,j)
+ end do
+ else if (k.eq.4) then
+ do j = 1,n_integ_pts + 1
+ integ(j) = Vn4_integ_save(i,pid,j)
+ xfunc(j) = Vn4_xfunc_save(i,pid,j)
+ end do
+ end if
+ Call LAGRNG(ran(irand),integ,Vn_event(i,k,pid),xfunc,
+ 1 n_integ_pts+1,1,5,n_integ_pts+1,1)
+ end if ! (For n_stdev_Vn*Vn_stdev=0, already have Vn_event)
+ end do
+ end do
+ end if
+
+ end do ! End Particle ID setup Loop
+
+ event_abort = 1
+
+ if(model_type .le. 4) then
+CCC Check validity of Temperature and sigma parameters:
+ do pid = 1,n_pid_type
+ if(Temp_event(pid) .le. Temp_abort) event_abort = -86
+ if(sigma_event(pid).le.sigma_abort) event_abort = -86
+ end do
+ end if
+
+ if(event_abort .eq. 1) then
+
+ total_mult = 0
+ do pid = 1,n_pid_type
+ total_mult = total_mult + mult_event(pid)
+ end do
+ n_vertices = 0
+ status_abort = 1
+ do pid = 1,n_pid_type
+CCC Load Anisotropic flow parameters for this event# and PID type
+CCC into temporary array;
+ do i = 1,nflowterms
+ do k = 1,4
+ Vn_event_tmp(i,k) = Vn_event(i,k,pid)
+ end do
+ end do
+CCC For the specified Geant PID, multiplicity, model type, temperature,
+CCC rapidity width (sigma), and expansion velocity parameter, generate
+CCC random track list:
+
+ Call track_gen(gpid(pid),mult_event(pid),model_type,
+ 1 Temp_event(pid),sigma_event(pid),expvel_event(pid),
+ 2 reac_plane_cntrl,PSIr_event,Vn_event_tmp,
+ 3 pt_cut_min,pt_cut_max,eta_cut_min,eta_cut_max,
+ 4 y_cut_min,y_cut_max,
+ 5 phi_cut_min,phi_cut_max,n_scan_pts,irand,rad,pid,status,
+ 6 n_integ_pts)
+ if(status .eq. -86) status_abort = -86
+ end do
+ end if
+
+ if(event_abort.eq.1 .and. status_abort.eq.1) then
+CCC Event Header Output:
+C--> write(7,230) ievent, total_mult
+C--> write(7,2301) PSIr_event
+C--> write(7,2302) MultFac_event
+C--> write(7,213) (gpid(pid),pid = 1,n_pid_type)
+C--> write(7,231) (mult_event(pid),pid = 1,n_pid_type)
+C--> if(model_type .le. 4) then
+C--> write(7,232) (Temp_event(pid),pid = 1,n_pid_type)
+C--> write(7,233) (sigma_event(pid),pid = 1,n_pid_type)
+C--> write(7,234) (expvel_event(pid),pid = 1,n_pid_type)
+C--> end if
+
+230 Format(/'*** Next Event: No. ',I7,'; Total # tracks = ',I10)
+2301 Format('* Reaction plane angle = ',G12.5,' (deg)')
+2302 Format('* Multiplicity Scaling Factor = ',G12.5)
+231 Format('* Multiplicity: ',10I7)
+232 Format('* Temperature: ',10F7.4)
+233 Format('* Rapidity Dist. sigma: ',10F7.4)
+234 Format('* Expansion Velocity: ',10F7.4)
+
+CCC Print out flow parameters for this event:
+C--> do pid = 1,n_pid_type
+C--> do i = 1,nflowterms
+C--> write(7,2341) pid,i,(Vn_event(i,k,pid),k=1,4)
+C--> end do
+C--> end do
+2341 Format('* Anisotropy Vn (pid=',I2,',n=',I1,'):',4G12.5)
+
+C--> write(7,235) ievent,total_mult,n_vertices
+235 Format('EVENT:',3x,3(1x,I6))
+C--> write(7,236)
+C--> write(7,237)
+236 Format('*** Tracks for this event ***')
+237 Format('* Geant PID # px py pz')
+CCC End Event Header Output:
+
+CCC Output track kinematics for ievent and pid:
+ track_counter = 0
+ start_v = 0
+ stop_v = 0
+ do pid = 1,n_pid_type
+ do i = 1,mult_event(pid)
+ track_counter = track_counter + 1
+C--> write(7,240) gpid(pid),(pout(pid,j,i),j=1,3),
+C--> 1 track_counter,start_v,stop_v,gpid(pid)
+C-OUT masstemp = pout(pid,4,i)
+C-OUT pxtemp = pout(pid,1,i)
+C-OUT pytemp = pout(pid,2,i)
+C-OUT pztemp = pout(pid,3,i)
+C-OUT Call kinematics2(masstemp,pxtemp,pytemp,pztemp,pttemp,
+C-OUT 1 etatemp,ytemp,phitemp)
+C-OUT write(9,270) gpid(pid),masstemp,pttemp,etatemp,ytemp,phitemp
+C-OUT270 Format(1x,I4,5E15.6)
+ masstemp = pout(pid,4,i)
+ pxtemp = pout(pid,1,i)
+ pytemp = pout(pid,2,i)
+ pztemp = pout(pid,3,i)
+ Call kinematics2(masstemp,pxtemp,pytemp,pztemp,pttemp,
+ 1 etatemp,ytemp,phitemp)
+C--> write(9,270) gpid(pid),masstemp,pttemp,etatemp,ytemp,phitemp
+270 Format(1x,I4,5E15.6)
+C-OUT Compute the cos(n*phi) Fourier component projections of the
+C-OUT azimuthal distributions for each PID type:
+ total_mult_inc(pid) = total_mult_inc(pid) + 1
+ jodd = 1
+ do j = 1,nflowterms
+ if(jodd.eq.1) then
+ if(ytemp.ge.0.0) then
+ cosinefac(j) =
+ 1 cos(float(j)*(phitemp-PSIr_event)/rad)
+ else if(ytemp.lt.0.0) then
+ cosinefac(j) =
+ 1 -cos(float(j)*(phitemp-PSIr_event)/rad)
+ end if
+ else if(jodd.eq.(-1)) then
+ cosinefac(j) =
+ 1 cos(float(j)*(phitemp-PSIr_event)/rad)
+ end if
+ jodd = -jodd
+ Vn_sum(j,pid) = Vn_sum(j,pid) + cosinefac(j)
+ end do
+C--> write(10,260) gpid(pid),(cosinefac(j),j=1,nflowterms)
+260 Format(1x,I3,6E12.4)
+C-OUT
+
+ end do
+ end do
+240 Format('TRACK:',1x,I6,3(1x,G12.5),4(1x,I6))
+CCC End track kinematics output.
+
+ else
+C--> write(7,250) ievent, event_abort, status_abort
+C--> write(8,250) ievent, event_abort, status_abort
+250 Format(/'*** Event No. ',I7,' ABORTED: event_abort,',
+ 1 'status_abort = ',2I7)
+ end if
+
+ END DO ! End Event Loop
+
+CCC Output Anisotropic flow check sums to terminal:
+ do pid = 1,n_pid_type
+C--> write(6,*) 'Total incl # part type (',gpid(pid),
+C--> 1 ') = ',total_mult_inc(pid)
+ do j = 1,nflowterms
+ Vn_sum(j,pid) = Vn_sum(j,pid)/total_mult_inc(pid)
+C--> write(6,*) 'Flow term #',j,': Check sum = ',
+C--> 1 Vn_sum(j,pid)
+ end do
+ end do
+
+CCC Output File Termination:
+ ievent = -999
+ total_mult = 0
+ n_vertices = 0
+C--> write(7,241)
+C--> write(7,235) ievent,total_mult,n_vertices
+241 Format(/'*** End of File ***')
+
+ Close(Unit=4)
+ Close(Unit=7)
+ Close(Unit=8)
+C-OUT Close(Unit=9)
+C-OUT Close(Unit=10)
+ Close(Unit=9)
+ Close(Unit=10)
+ RETURN
+ END
+
+ Subroutine track_gen(gpid,N,model_type,T,sigma,expvel,
+ 1 reac_plane_cntrl,PSIr,Vn,
+ 2 pt_cut_min,pt_cut_max,eta_cut_min,eta_cut_max,
+ 3 y_cut_min,y_cut_max,
+ 4 phi_cut_min,phi_cut_max,n_scan_pts,irand,rad,pid,status,
+ 5 n_integ_pts)
+ implicit none
+ Include 'common_facfac.inc'
+ Include 'Parameter_values.inc'
+ Common/track/ pout(npid,4,factorial_max)
+ real*4 pout
+ Common/Geant/geant_mass(200),geant_charge(200),
+ 1 geant_lifetime(200),geant_width(200)
+ real*4 geant_mass,geant_charge,geant_lifetime,geant_width
+
+ real*4 T,sigma,expvel,pt_cut_min,pt_cut_max,eta_cut_min
+ real*4 eta_cut_max,phi_cut_min,phi_cut_max,mass
+ real*4 dpt,deta,facmax,pt,eta,y,rapidity,dNdpty,dNdp
+ real*4 pt_trial,eta_trial,y_trial,ran,rad,phi
+ real*4 y_cut_min,y_cut_max,pseudorapidity
+ real*4 PSIr,Vn(nflowterms,4),dphi,facmax_phi,dNdphi
+ real*4 Vnptmax,Vnymax,Vn_pt_y,Vntmp(nflowterms)
+ real*4 masstmp,Mass_function
+
+ integer gpid,N,model_type,npt,neta,n_scan_pts,ipt,ieta
+ integer Track_count,irand,i,j,pid,status
+ integer reac_plane_cntrl,iphi
+ integer n_integ_pts
+
+ do i = 1,factorial_max
+ do j = 1,4
+ pout(pid,j,i) = 0.0
+ end do
+ end do
+
+ mass = geant_mass(gpid)
+ npt = n_scan_pts
+ neta = n_scan_pts
+
+CCC Determine maximum value of model distribution in (pt,eta) range:
+
+ dpt = (pt_cut_max - pt_cut_min)/float(npt - 1)
+ deta = (eta_cut_max - eta_cut_min)/float(neta - 1)
+ facmax = 0.0
+ do ipt = 1,npt
+ pt = pt_cut_min + dpt*float(ipt - 1)
+ do ieta = 1,neta
+ eta = eta_cut_min + deta*float(ieta - 1)
+ y = rapidity(mass,pt,eta)
+ dNdp = dNdpty(1.0,pt,eta,y,mass,T,sigma,expvel,
+ 1 model_type,1,pid)
+ if(dNdp .gt. facmax) facmax = dNdp
+ end do
+ end do
+
+CCC If dNdp always underflows exp() range, then facmax will stay
+CCC equal to 0.0, thus causing a divide by 0.0 error below.
+CCC Check for this and Return if this is the case. This event will
+CCC be aborted in this instance.
+
+ if(facmax .eq. 0.0) then
+ status = -86
+ Return
+ else
+ status = 1
+ end if
+
+CCC Determine the maximum values of the azimuthal model distribution
+CCC in phi, for case with reaction plane dependence and anisotropic flow
+CCC Find the absolute maximum possible value given the pt and y dependences
+CCC and assuming all Fourier components add with maximum coherence.
+
+ if(reac_plane_cntrl .gt. 1) then
+ facmax_phi = 1.0
+ do i = 1,nflowterms
+ if(i.eq.(2*(i/2))) then ! Select even Fourier components:
+ Vnptmax = max(
+ 1 abs(Vn(i,1)+Vn(i,4)*pt_cut_min
+ 3 +Vn(i,2)*pt_cut_min*pt_cut_min),
+ 2 abs(Vn(i,1)+Vn(i,4)*pt_cut_max
+ 4 +Vn(i,2)*pt_cut_max*pt_cut_max))
+ Vnymax = max(
+ 1 exp(-Vn(i,3)*y_cut_min*y_cut_min),
+ 2 exp(-Vn(i,3)*y_cut_max*y_cut_max))
+ if(y_cut_min.lt.0.0 .and. y_cut_max.gt.0.0) then
+ Vnymax = max(Vnymax,1.0)
+ end if
+ else ! Select odd Fourier components
+ Vnptmax = max(
+ 1 abs(Vn(i,1)+Vn(i,2)*pt_cut_min),
+ 2 abs(Vn(i,1)+Vn(i,2)*pt_cut_max))
+ Vnymax = max(
+ 1 abs(Vn(i,3)+Vn(i,4)*abs(y_cut_min**3)),
+ 2 abs(Vn(i,3)+Vn(i,4)*abs(y_cut_max**3)))
+ if(y_cut_min.lt.0.0 .and. y_cut_max.gt.0.0) then
+ Vnymax = max(Vnymax,abs(Vn(i,3)))
+ end if
+ end if
+ facmax_phi = facmax_phi + 2.0*Vnptmax*Vnymax
+ end do
+CCC Check facmax_phi; if 0, then event will be aborted.
+ if(facmax_phi.eq.0.0) then
+ status = -86
+ Return
+ else
+ status = 1
+ end if
+ end if
+
+CCC Start Random Track Selection:
+
+ Track_count = 0
+
+100 pt_trial = ran(irand)*(pt_cut_max - pt_cut_min)+pt_cut_min
+ if(model_type.ge.1 .and. model_type.le.4) then
+ y_trial = ran(irand)*(y_cut_max - y_cut_min)+y_cut_min
+ eta_trial = pseudorapidity(mass,pt_trial,y_trial)
+ if(eta_trial.lt.eta_cut_min .or. eta_trial.gt.eta_cut_max)
+ 1 go to 100
+ else if (model_type.eq.5 .or. model_type.eq.6) then
+ eta_trial=ran(irand)*(eta_cut_max-eta_cut_min)+eta_cut_min
+ y_trial = rapidity(mass,pt_trial,eta_trial)
+ end if
+ dNdp = dNdpty(1.0,pt_trial,eta_trial,y_trial,mass,T,sigma,
+ 1 expvel,model_type,1,pid)/facmax
+ if(ran(irand) .le. dNdp) then
+ Track_count = Track_count + 1
+
+CCC Determine phi angle:
+ if(reac_plane_cntrl .eq. 1) then
+ phi = (ran(irand)*(phi_cut_max - phi_cut_min) +
+ 1 phi_cut_min)/rad
+ else if(reac_plane_cntrl .gt. 1) then
+ do i = 1,nflowterms
+ Vntmp(i) = Vn_pt_y(i,Vn(i,1),Vn(i,2),Vn(i,3),Vn(i,4),
+ 1 pt_trial,y_trial)
+ end do
+101 phi = ran(irand)*(phi_cut_max - phi_cut_min)+phi_cut_min
+ dNdphi = 1.0
+ do i = 1,nflowterms
+ dNdphi=dNdphi+2.0*Vntmp(i)*cos(float(i)*(phi-PSIr)/rad)
+ end do
+ if(ran(irand).gt.(dNdphi/facmax_phi)) then
+ go to 101
+ else
+ phi = phi/rad
+ end if
+ end if
+
+ masstmp = Mass_function(gpid,pid,n_integ_pts,irand)
+ Call kinematics(masstmp,pt_trial,y_trial,phi,Track_count,pid)
+ if(Track_count .lt. N) then
+ go to 100
+ else if(Track_count .ge. N) then
+ Return
+ end if
+
+ else
+ go to 100
+ end if
+
+ END
+
+ Real*4 Function rapidity(m,pt,eta)
+ implicit none
+ real*4 m,pt,eta,theta,pz,E,pi,expR
+
+ pi = 3.141592654
+ theta = 2.0*ATAN(expR(-eta))
+ if(abs(theta - pi/2.0) .lt. 0.000001) then
+ pz = 0.0
+ else
+ pz = pt/tan(theta)
+ end if
+ E = sqrt(pt*pt + pz*pz + m*m)
+ rapidity = 0.5*log((E+pz)/(E-pz))
+ Return
+ END
+
+ Real*4 Function pseudorapidity(m,pt,y)
+ implicit none
+
+CCC This Function computes the pseudorapidty for a given mass, pt, y:
+
+ real*4 m,pt,y,mt,coshy,E,pzmag,pz,pmag,expR
+
+ if(y.eq.0.0) then
+ pseudorapidity = 0.0
+ Return
+ end if
+ mt = sqrt(m*m + pt*pt)
+ coshy = 0.5*(expR(y) + expR(-y))
+ E = mt*coshy
+ pzmag = sqrt(abs(E*E - mt*mt))
+ if(y.eq.0.0) then
+ pz = 0.0
+ else
+ pz = (y/abs(y))*pzmag
+ end if
+ pmag = sqrt(pt*pt + pz*pz)
+ if(pt.ne.0.0) then
+ pseudorapidity = 0.5*log((pmag+pz)/(pmag-pz))
+ else if (pt.eq.0.0) then
+ pseudorapidity = 86.0
+C--> write(8,10)
+10 Format(10x,'Function pseudorapidity called with pt=0')
+ end if
+ Return
+ END
+
+ Subroutine kinematics(m,pt,y,phi,index,pid)
+ implicit none
+
+CCC This SUBR takes the input particle mass (m), pt, y and phi and
+CCC computes E, px, py, pz and loads the momenta into the index-th
+CCC row of array pout(,,) in Common/track/ .
+
+ integer index,pid
+ Include 'common_facfac.inc'
+ Include 'Parameter_values.inc'
+ Common/track/ pout(npid,4,factorial_max)
+ real*4 pout
+
+ real*4 m,pt,y,phi,mt,coshy,E,pzmag,px,py,pz,expR
+
+ mt = sqrt(m*m + pt*pt)
+ coshy = 0.5*(expR(y) + expR(-y))
+ E = mt*coshy
+ pzmag = sqrt(abs(E*E - mt*mt))
+ if(y.eq.0.0) then
+ pz = 0.0
+ else
+ pz = (y/abs(y))*pzmag
+ end if
+ px = pt*cos(phi)
+ py = pt*sin(phi)
+
+ pout(pid,1,index) = px
+ pout(pid,2,index) = py
+ pout(pid,3,index) = pz
+ pout(pid,4,index) = m
+
+ Return
+ END
+
+ Subroutine kinematics2(m,px,py,pz,pt,eta,y,phi)
+ implicit none
+
+CCC This SUBR takes the input particle mass (m), px,py,pz and
+CCC computes pt,eta,y,phi:
+
+ real*4 m,px,py,pz,pt,eta,y,phi,pi,E,E0
+
+ pi = 3.141592654
+ pt = sqrt(px*px + py*py)
+ E = sqrt(m*m + pz*pz + pt*pt)
+ y = 0.5*log((E + pz)/(E - pz))
+ E0 = sqrt(pz*pz + pt*pt)
+ if(pt.eq.0.0) then
+ eta = -86.0
+ else
+ eta = 0.5*log((E0 + pz)/(E0 - pz))
+ end if
+ if(px.eq.0.0 .and. py.eq.0.0) then
+ phi = 0.0
+ else
+ phi = atan2(py,px)
+ end if
+ phi = (180.0/pi)*phi
+ if(phi .lt. 0.0) phi = phi + 360.0
+ Return
+ END
+
+ Real*4 Function dNdpty(A,pt,eta,y,m,T,sigma,vel,control,ktl,pid)
+ implicit none
+
+ real*4 A,pt,eta,y,m,T,sigma,vel
+ real*4 pi,mt,coshy,E,ptot,FAC,gamma,yp,sinhyp,coshyp
+ real*4 FAC2,FAC3,expR
+ integer control, ktl,pid,pt_index,eta_index,index_locate
+ Include 'Parameter_values.inc'
+ Include 'common_dist_bin.inc'
+
+CCC Calculates dN/dp^3 using several models:
+CCC
+CCC control = 1, Humanic factorized model
+CCC = 2, Pratt non-expanding spherical thermal source
+CCC = 3, Bertsch non-expanding spherical thermal source
+CCC = 4, Pratt spherical expanding thermally equilibrated
+CCC source.
+CCC = 5, Factorized pt, eta bin-by-bin distribution.
+CCC = 6, Full 2D pt, eta bin-by-bin distribution.
+CCC
+CCC ktl = 0, to return value of dN/dp^3
+CCC ktl = 1, to return value of dN/dpt*dy
+
+ pi = 3.141592654
+ mt = sqrt(pt*pt + m*m)
+ coshy = 0.5*(expR(y) + expR(-y))
+ E = mt*coshy
+ ptot = sqrt(E*E - m*m)
+ if(ktl .eq. 0) then
+ FAC = 2.0*pi*pt*E
+ else if(ktl .eq. 1) then
+ FAC = 1.0
+ end if
+
+ if(control .eq. 1) then
+ dNdpty = A*pt*expR(-mt/T)*expR(-y*y/(2.0*sigma*sigma))
+ dNdpty = dNdpty/FAC
+
+ else if(control .eq. 2) then
+ dNdpty = A*pt*E*expR(-E/T)
+ dNdpty = dNdpty/FAC
+
+ else if(control .eq. 3) then
+ dNdpty = A*pt*E/(expR(E/T) - 1.0)
+ dNdpty = dNdpty/FAC
+
+ else if(control .eq. 4) then
+ gamma = 1.0/sqrt(1.0 - vel*vel)
+ yp = gamma*vel*ptot/T
+ sinhyp = 0.5*(expR(yp) - expR(-yp))
+ coshyp = 0.5*(expR(yp) + expR(-yp))
+ if(yp .ne. 0.0) then
+ FAC2 = sinhyp/yp
+ else
+ FAC2 = 1.0
+ end if
+ FAC3 = FAC2+ (T/(gamma*E))*(FAC2 - coshyp)
+ dNdpty = A*pt*E*expR(-gamma*E/T)*FAC3
+ dNdpty = dNdpty/FAC
+
+ else if(control .eq. 5) then
+ pt_index = index_locate(pid,pt,1)
+ eta_index = index_locate(pid,eta,2)
+ dNdpty = A*pt_bin(pid,pt_index)*eta_bin(pid,eta_index)
+ dNdpty = dNdpty/FAC
+
+ else if(control .eq. 6) then
+ pt_index = index_locate(pid,pt,1)
+ eta_index = index_locate(pid,eta,2)
+ dNdpty = A*pt_eta_bin(pid,pt_index,eta_index)
+ dNdpty = dNdpty/FAC
+
+ else
+ dNdpty = -86.0
+
+ end if
+
+ return
+ END
+
+ Integer Function index_locate(pid,arg,kind)
+ implicit none
+
+ Include 'Parameter_values.inc'
+ Include 'common_dist_bin.inc'
+
+ integer pid,kind,ibin
+ real*4 arg
+
+CCC This Function locates the pt or eta bin number corresponding to the
+CCC input bin mesh, the requested value of pt or eta, for the current
+CCC value of particle type.
+CCC
+CCC If kind = 1, then pt bin number is located.
+CCC If kind = 2, then eta bin number is located.
+
+ if(kind .eq. 1) then
+ do ibin = 1,n_pt_bins(pid)
+ if(arg.le.pt_bin_mesh(pid,ibin)) then
+ index_locate = ibin
+ Return
+ end if
+ end do
+ index_locate = n_pt_bins(pid)
+C--> write(8,10) pid,arg
+10 Format(//10x,'In Function index_locate, for pid = ',I5,
+ 1 'pt =',E15.6,' is out of range - use last bin #')
+ Return
+
+ else if(kind .eq. 2) then
+
+ do ibin = 1,n_eta_bins(pid)
+ if(arg.le.eta_bin_mesh(pid,ibin)) then
+ index_locate = ibin
+ Return
+ end if
+ end do
+ index_locate = n_eta_bins(pid)
+C--> write(8,11) pid,arg
+11 Format(//10x,'In Function index_locate, for pid = ',I5,
+ 1 'eta =',E15.6,' is out of range - use last bin #')
+ Return
+
+ end if
+ END
+
+ Real*4 Function expR(x)
+ implicit none
+ real*4 x
+ if(x .gt. 69.0) then
+C--> write(8,10) x
+ STOP
+ else if(x .lt. -69.0) then
+ expR = 0.0
+ else
+ expR = exp(x)
+ end if
+10 Format(///10x,'Func. expR(x) called with x = ',E15.6,
+ 1 ', gt 69.0 - STOP')
+ Return
+ END
+
+ SUBROUTINE LAGRNG (X,ARG,Y,VAL,NDIM,NFS,NPTS,MAXARG,MAXFS)
+ IMPLICIT REAL*4(A-H,O-Z)
+C
+C LAGRANGE INTERPOLATION,UNEQUALLY SPACED POINTS
+C ROUTINE OBTAINED FROM R. LANDAU, UNIV. OF OREGON.
+C ARG=VECTOR OF INDEPENDENT VARIABLE CONTAINING MAXARG VALUES.
+C VAL=MATRIX OF FUNCTION VALUES CORRESPONDING TO ARG. (MAXFS
+C FUNCTIONS AT MAXARG VALUES.)
+C X =VALUE OF INDEP. VARIABLE FOR WHICH INTERPOLATION IS DESIRED.
+C Y =VECTOR OF MAXFS FUNCTION VALUES RESULTING FROM SIMUL. INTERP.
+C NDIM=NUMBER OF ARG VALUES TO BE USED. (NDIM.LE.MAXARG)
+C NFS=NUMBER OF FUNCTIONS SIMUL. INTERP (NFS.LE.MAXFS)
+C NPTS=NUMBER OF POINTS USED IN INTERPOLATION. (NPTS=2,3,4,5,6)
+C
+ DIMENSION ARG(MAXARG), VAL(MAXFS,MAXARG), Y(MAXFS)
+C
+C -----FIND X0, THE CLOSEST POINT TO X.
+C
+ NI=1
+ NF=NDIM
+ 10 IF ((X.LE.ARG(NI)).OR.(X.GE.ARG(NF))) GO TO 30
+ IF ((NF-NI+1).EQ.2) GO TO 70
+ NMID=(NF+NI)/2
+ IF (X.GT.ARG(NMID)) GO TO 20
+ NF=NMID
+ GO TO 10
+ 20 NI=NMID
+ GO TO 10
+C
+C ------ X IS ONE OF THE TABLULATED VALUES.
+C
+ 30 IF (X.LE.ARG(NI)) GO TO 60
+ NN=NF
+ 40 NUSED=0
+ DO 50 N=1,NFS
+ 50 Y(N)=VAL(N,NN)
+ RETURN
+ 60 NN=NI
+ GO TO 40
+C
+C ------- 2 PTS LEFT, CHOOSE SMALLER ONE.
+C
+ 70 N0=NI
+ NN=NPTS-2
+ GO TO (110,100,90,80), NN
+ 80 CONTINUE
+ IF (((N0+3).GT.NDIM).OR.((N0-2).LT.1)) GO TO 90
+ NUSED=6
+ GO TO 130
+ 90 CONTINUE
+ IF ((N0+2).GT.NDIM) GO TO 110
+ IF ((N0-2).LT.1) GO TO 100
+ NUSED=5
+ GO TO 130
+ 100 CONTINUE
+ IF (((N0+2).GT.NDIM).OR.((N0-1).LT.1)) GO TO 110
+ NUSED=4
+ GO TO 130
+ 110 IF ((N0+1).LT.NDIM) GO TO 120
+C
+C ------N0=NDIM, SPECIAL CASE.
+C
+ NN=NDIM
+ GO TO 40
+ 120 NUSED=3
+ IF ((N0-1).LT.1) NUSED=2
+ 130 CONTINUE
+C
+C ------AT LEAST 2 PTS LEFT.
+C
+ Y0=X-ARG(N0)
+ Y1=X-ARG(N0+1)
+ Y01=Y1-Y0
+ C0=Y1/Y01
+ C1=-Y0/Y01
+ IF (NUSED.EQ.2) GO TO 140
+C
+C ------AT LEAST 3 PTS.
+C
+ YM1=X-ARG(N0-1)
+ Y0M1=YM1-Y0
+ YM11=Y1-YM1
+ CM1=-Y0*Y1/Y0M1/YM11
+ C0=C0*YM1/Y0M1
+ C1=-C1*YM1/YM11
+ IF (NUSED.EQ.3) GO TO 160
+C
+C ------AT LEAST 4 PTS
+C
+ Y2=X-ARG(N0+2)
+ YM12=Y2-YM1
+ Y02=Y2-Y0
+ Y12=Y2-Y1
+ CM1=CM1*Y2/YM12
+ C0=C0*Y2/Y02
+ C1=C1*Y2/Y12
+ C2=-YM1*Y0*Y1/YM12/Y02/Y12
+ IF (NUSED.EQ.4) GO TO 180
+C
+C ------AT LEAST 5 PTS.
+C
+ YM2=X-ARG(N0-2)
+ YM2M1=YM1-YM2
+ YM20=Y0-YM2
+ YM21=Y1-YM2
+ YM22=Y2-YM2
+ CM2=YM1*Y0*Y1*Y2/YM2M1/YM20/YM21/YM22
+ CM1=-CM1*YM2/YM2M1
+ C0=-C0*YM2/YM20
+ C1=-C1*YM2/YM21
+ C2=-C2*YM2/YM22
+ IF (NUSED.EQ.5) GO TO 200
+C
+C ------AT LEAST 6 PTS.
+C
+ Y3=X-ARG(N0+3)
+ YM23=Y3-YM2
+ YM13=Y3-YM1
+ Y03=Y3-Y0
+ Y13=Y3-Y1
+ Y23=Y3-Y2
+ CM2=CM2*Y3/YM23
+ CM1=CM1*Y3/YM13
+ C0=C0*Y3/Y03
+ C1=C1*Y3/Y13
+ C2=C2*Y3/Y23
+ C3=YM2*YM1*Y0*Y1*Y2/YM23/YM13/Y03/Y13/Y23
+ GO TO 220
+ 140 CONTINUE
+ DO 150 N=1,NFS
+ 150 Y(N)=C0*VAL(N,N0)+C1*VAL(N,N0+1)
+ GO TO 240
+ 160 CONTINUE
+ DO 170 N=1,NFS
+ 170 Y(N)=CM1*VAL(N,N0-1)+C0*VAL(N,N0)+C1*VAL(N,N0+1)
+ GO TO 240
+ 180 CONTINUE
+ DO 190 N=1,NFS
+ 190 Y(N)=CM1*VAL(N,N0-1)+C0*VAL(N,N0)+C1*VAL(N,N0+1)+C2*VAL(N,N0+2)
+ GO TO 240
+ 200 CONTINUE
+ DO 210 N=1,NFS
+ 210 Y(N)=CM2*VAL(N,N0-2)+CM1*VAL(N,N0-1)+C0*VAL(N,N0)+C1*VAL(N,N0+1)+C
+ 12*VAL(N,N0+2)
+ GO TO 240
+ 220 CONTINUE
+ DO 230 N=1,NFS
+ 230 Y(N)=CM2*VAL(N,N0-2)+CM1*VAL(N,N0-1)+C0*VAL(N,N0)+C1*VAL(N,N0+1)+C
+ 12*VAL(N,N0+2)+C3*VAL(N,N0+3)
+ 240 RETURN
+C
+ END
+
+ Subroutine FACTORIAL
+ implicit none
+
+ integer n
+ Include 'common_facfac.inc'
+ real*4 FN
+
+CCC Computes the natural log of n! for n = 0,1,2...(factorial_max -1)
+CCC and puts the result in array FACLOG().
+C
+CCC FACLOG(1) = log(0!) = 0.0
+CCC FACLOG(2) = log(1!) = 0.0
+CCC FACLOG(n+1) = log(n!)
+
+ FACLOG(1) = 0.0
+ FACLOG(2) = 0.0
+ FN = 1.0
+ do n = 3,factorial_max
+ FN = FN + 1.0
+ FACLOG(n) = FACLOG(n-1) + log(FN)
+ end do
+ Return
+ END
+
+ Subroutine MinMax(mean,stdev,n_stdev,min,max)
+ implicit none
+
+CCC Computes range of integration for random number selections
+
+ real*4 mean,stdev,n_stdev,min,max
+
+ min = mean - n_stdev*stdev
+ if(min .lt. 0.0) min = 0.0
+ max = mean + n_stdev*stdev
+ Return
+ END
+
+ Subroutine Poisson(min,max,mean,nsteps,integ,xfunc,ndim)
+ implicit none
+
+CCC Computes Poisson distribution from n = min to max;
+CCC Integrates this distribution and records result at each step in
+CCC array integ();
+CCC Records the coordinates in array xfunc().
+
+ integer min,max,mean,nsteps,ndim,i,n
+ Include 'Parameter_values.inc'
+ real*4 mean_real,mean_real_log,expR
+ real*4 integ(ndim)
+ real*4 xfunc(ndim)
+ real*4 Poisson_dist(n_mult_max_steps)
+ Include 'common_facfac.inc'
+
+CCC Initialize arrays to zero:
+
+ do i = 1,ndim
+ integ(i) = 0.0
+ xfunc(i) = 0.0
+ Poisson_dist(i) = 0.0
+ end do
+
+ mean_real = float(mean)
+ mean_real_log = log(mean_real)
+
+CCC Compute Poisson distribution from n = min to max:
+ do i = 1,nsteps
+ n = (i - 1) + min
+ Poisson_dist(i) = expR(-mean_real + float(n)*mean_real_log
+ 1 - FACLOG(n+1))
+ end do
+
+CCC Integrate the Poisson distribution:
+ integ(1) = 0.0
+ do i = 2,nsteps
+ integ(i) = 0.5*(Poisson_dist(i-1) + Poisson_dist(i)) + integ(i-1)
+ end do
+
+CCC Normalize the integral to unity:
+ do i = 1,nsteps
+ integ(i) = integ(i)/integ(nsteps)
+ end do
+
+CCC Fill xfunc array:
+ do i = 1,nsteps
+ xfunc(i) = float(i - 1 + min)
+ end do
+
+CCC Extend integ() and xfunc() by one additional mesh point past the
+CCC end point in order to avoid a bug in the Lagrange interpolation
+CCC subroutine that gives erroneous interpolation results within the
+CCC the last mesh bin.
+
+ integ(nsteps + 1) = integ(nsteps) + 0.01
+ xfunc(nsteps + 1) = xfunc(nsteps)
+
+ Return
+ END
+
+ Subroutine Gaussian(min,max,mean,stdev,npts,integ,xfunc,ndim)
+ implicit none
+
+CCC Compute Gaussian distribution from x = min to max at npts;
+CCC Integrate this distribution and record result at each mesh in
+CCC array integ();
+CCC Record the coordinates in array xfunc().
+
+ Include 'Parameter_values.inc'
+ integer npts,ndim,i
+ real*4 min,max,mean,stdev,integ(ndim),xfunc(ndim)
+ real*4 dm,x,Gauss_dist(nmax_integ),FAC1,FAC2,pi,expR
+
+CCC Initialize arrays to zero:
+ do i = 1,ndim
+ integ(i) = 0.0
+ xfunc(i) = 0.0
+ Gauss_dist(i) = 0.0
+ end do
+
+ pi = 3.141592654
+ FAC1 = 1.0/(sqrt(2.0*pi)*stdev)
+ FAC2 = 2.0*stdev*stdev
+ dm = (max - min)/float(npts - 1)
+
+CCC Compute normalized Gaussian distribution:
+ do i = 1,npts
+ x = min + dm*float(i-1)
+ xfunc(i) = x
+ Gauss_dist(i) = FAC1*expR(-((x - mean)**2)/FAC2)
+ end do
+
+CCC Integrate Gaussian distribution over specified range
+ integ(1) = 0.0
+ do i = 2,npts
+ integ(i) = 0.5*(Gauss_dist(i-1) + Gauss_dist(i))*dm + integ(i-1)
+ end do
+
+CCC Normalize integral to unity:
+ do i = 1,npts
+ integ(i) = integ(i)/integ(npts)
+ end do
+
+CCC Extend integ() and xfunc() by one mesh point to avoid Lagrange
+CCC interpolation subroutine bug:
+ integ(npts + 1) = integ(npts) + 0.01
+ xfunc(npts + 1) = xfunc(npts)
+
+ Return
+ END
+
+ Subroutine Particle_prop
+ implicit none
+
+ Common/Geant/geant_mass(200),geant_charge(200),
+ 1 geant_lifetime(200),geant_width(200)
+
+CCC Local Variable Type Declarations:
+ integer i
+ real*4 geant_mass,geant_charge,geant_lifetime,geant_width
+ real*4 hbar
+ Parameter(hbar = 6.5821220E-25) ! hbar in GeV*seconds
+
+CCC Fill arrays geant_mass(),geant_charge(),geant_lifetime(),
+CCC geant_width() with particle mass in GeV, charge in units of
+CCC |e|, mean lifetime in seconds, and width in GeV, where
+CCC width*lifetime = hbar. For lifetimes listed with values of
+CCC 1.0E+30 (i.e. stable particles) the width is set to 0.000.
+CCC Row # = Geant Particle ID # code (see Geant Manual 3.10,
+CCC User's Guide, CONS 300-1 and 300-2). Note, rows 151 and higher
+CCC are used for resonances. These assume the masses and widths
+CCC specific to the models used to represent the invariant mass
+CCC distributions and therefore may differ slightly from the Particle
+CCC Data Group's listing.
+
+CCC Initialize arrays to zero:
+ do i = 1,200
+ geant_mass(i) = 0.0
+ geant_charge(i) = 0.0
+ geant_lifetime(i) = 0.0
+ geant_width(i) = 0.0
+ end do
+
+CCC Set Particle Masses in GeV:
+ geant_mass(1) = 0.0 ! Gamma
+ geant_mass(2) = 0.000511 ! Positron
+ geant_mass(3) = 0.000511 ! Electron
+ geant_mass(4) = 0.0 ! Neutrino
+ geant_mass(5) = 0.105659 ! Muon +
+ geant_mass(6) = 0.105659 ! Muon -
+ geant_mass(7) = 0.134693 ! Pion 0
+ geant_mass(8) = 0.139567 ! Pion +
+ geant_mass(9) = 0.139567 ! Pion -
+ geant_mass(10) = 0.49767 ! Kaon 0 Long
+ geant_mass(11) = 0.493667 ! Kaon +
+ geant_mass(12) = 0.493667 ! Kaon -
+ geant_mass(13) = 0.939573 ! Neutron
+ geant_mass(14) = 0.93828 ! Proton
+ geant_mass(15) = 0.93828 ! Antiproton
+ geant_mass(16) = 0.49767 ! Kaon 0 Short
+ geant_mass(17) = 0.5488 ! Eta
+ geant_mass(18) = 1.11560 ! Lambda
+ geant_mass(19) = 1.18936 ! Sigma +
+ geant_mass(20) = 1.19246 ! Sigma 0
+ geant_mass(21) = 1.19734 ! Sigma -
+ geant_mass(22) = 1.31490 ! Xi 0
+ geant_mass(23) = 1.32132 ! Xi -
+ geant_mass(24) = 1.67245 ! Omega
+ geant_mass(25) = 0.939573 ! Antineutron
+ geant_mass(26) = 1.11560 ! Antilambda
+ geant_mass(27) = 1.18936 ! Antisigma -
+ geant_mass(28) = 1.19246 ! Antisigma 0
+ geant_mass(29) = 1.19734 ! Antisigma +
+ geant_mass(30) = 1.3149 ! Antixi 0
+ geant_mass(31) = 1.32132 ! Antixi +
+ geant_mass(32) = 1.67245 ! Antiomega +
+ geant_mass(33) = 1.7842 ! Tau +
+ geant_mass(34) = 1.7842 ! Tau -
+ geant_mass(35) = 1.8694 ! D+
+ geant_mass(36) = 1.8694 ! D-
+ geant_mass(37) = 1.8647 ! D0
+ geant_mass(38) = 1.8647 ! Anti D0
+ geant_mass(39) = 1.9710 ! F+, now called Ds+
+ geant_mass(40) = 1.9710 ! F-, now called Ds-
+ geant_mass(41) = 2.2822 ! Lambda C+
+ geant_mass(42) = 80.8000 ! W+
+ geant_mass(43) = 80.8000 ! W-
+ geant_mass(44) = 92.9000 ! Z0
+ geant_mass(45) = 1.877 ! Deuteron
+ geant_mass(46) = 2.817 ! Tritium
+ geant_mass(47) = 3.755 ! Alpha
+ geant_mass(48) = 0.0 ! Geantino
+ geant_mass(49) = 2.80923 ! He3
+ geant_mass(50) = 0.0 ! Cherenkov
+ geant_mass(151) = 0.783 ! rho +
+ geant_mass(152) = 0.783 ! rho -
+ geant_mass(153) = 0.783 ! rho 0
+ geant_mass(154) = 0.782 ! omega 0
+ geant_mass(155) = 0.95750 ! eta'
+ geant_mass(156) = 1.0194 ! phi
+ geant_mass(157) = 3.09693 ! J/Psi
+ geant_mass(158) = 1.232 ! Delta -
+ geant_mass(159) = 1.232 ! Delta 0
+ geant_mass(160) = 1.232 ! Delta +
+ geant_mass(161) = 1.232 ! Delta ++
+ geant_mass(162) = 0.89183 ! K* +
+ geant_mass(163) = 0.89183 ! K* -
+ geant_mass(164) = 0.89610 ! K* 0
+
+CCC Set Particle Charge in |e|:
+ geant_charge( 1) = 0 ! Gamma
+ geant_charge( 2) = 1 ! Positron
+ geant_charge( 3) = -1 ! Electron
+ geant_charge( 4) = 0 ! Neutrino
+ geant_charge( 5) = 1 ! Muon+
+ geant_charge( 6) = -1 ! Muon-
+ geant_charge( 7) = 0 ! Pion0
+ geant_charge( 8) = 1 ! Pion+
+ geant_charge( 9) = -1 ! Pion-
+ geant_charge(10) = 0 ! Kaon 0 long
+ geant_charge(11) = 1 ! Kaon+
+ geant_charge(12) = -1 ! Kaon-
+ geant_charge(13) = 0 ! Neutron
+ geant_charge(14) = 1 ! Proton
+ geant_charge(15) = -1 ! Antiproton
+ geant_charge(16) = 0 ! Kaon 0 short
+ geant_charge(17) = 0 ! Eta
+ geant_charge(18) = 0 ! Lambda
+ geant_charge(19) = 1 ! Sigma+
+ geant_charge(20) = 0 ! Sigma0
+ geant_charge(21) = -1 ! Sigma-
+ geant_charge(22) = 0 ! Xi 0
+ geant_charge(23) = -1 ! Xi -
+ geant_charge(24) = -1 ! Omega
+ geant_charge(25) = 0 ! Antineutron
+ geant_charge(26) = 0 ! Antilambda
+ geant_charge(27) = -1 ! Anti-Sigma -
+ geant_charge(28) = 0 ! Anti-Sigma 0
+ geant_charge(29) = 1 ! Anti-Sigma +
+ geant_charge(30) = 0 ! AntiXi 0
+ geant_charge(31) = 1 ! AntiXi +
+ geant_charge(32) = 1 ! Anti-Omega +
+ geant_charge(33) = 1 ! Tau +
+ geant_charge(34) = -1 ! Tau -
+ geant_charge(35) = 1 ! D+
+ geant_charge(36) = -1 ! D-
+ geant_charge(37) = 0 ! D0
+ geant_charge(38) = 0 ! Anti D0
+ geant_charge(39) = 1 ! F+, now called Ds+
+ geant_charge(40) = -1 ! F-, now called Ds-
+ geant_charge(41) = 1 ! Lambda C+
+ geant_charge(42) = 1 ! W+
+ geant_charge(43) = -1 ! W-
+ geant_charge(44) = 0 ! Z0
+ geant_charge(45) = 1 ! Deuteron
+ geant_charge(46) = 1 ! Triton
+ geant_charge(47) = 2 ! Alpha
+ geant_charge(48) = 0 ! Geantino (Fake particle)
+ geant_charge(49) = 2 ! He3
+ geant_charge(50) = 0 ! Cerenkov (Fake particle)
+ geant_charge(151) = 1 ! rho+
+ geant_charge(152) = -1 ! rho-
+ geant_charge(153) = 0 ! rho 0
+ geant_charge(154) = 0 ! omega 0
+ geant_charge(155) = 0 ! eta'
+ geant_charge(156) = 0 ! phi
+ geant_charge(157) = 0 ! J/Psi
+ geant_charge(158) = -1 ! Delta -
+ geant_charge(159) = 0 ! Delta 0
+ geant_charge(160) = 1 ! Delta +
+ geant_charge(161) = 2 ! Delta ++
+ geant_charge(162) = 1 ! K* +
+ geant_charge(163) = -1 ! K* -
+ geant_charge(164) = 0 ! K* 0
+
+CCC Set Particle Lifetimes in seconds:
+ geant_lifetime( 1) = 1.0E+30 ! Gamma
+ geant_lifetime( 2) = 1.0E+30 ! Positron
+ geant_lifetime( 3) = 1.0E+30 ! Electron
+ geant_lifetime( 4) = 1.0E+30 ! Neutrino
+ geant_lifetime( 5) = 2.19703E-06 ! Muon+
+ geant_lifetime( 6) = 2.19703E-06 ! Muon-
+ geant_lifetime( 7) = 8.4E-17 ! Pion0
+ geant_lifetime( 8) = 2.603E-08 ! Pion+
+ geant_lifetime( 9) = 2.603E-08 ! Pion-
+ geant_lifetime(10) = 5.16E-08 ! Kaon 0 long
+ geant_lifetime(11) = 1.237E-08 ! Kaon+
+ geant_lifetime(12) = 1.237E-08 ! Kaon-
+ geant_lifetime(13) = 889.1 ! Neutron
+ geant_lifetime(14) = 1.0E+30 ! Proton
+ geant_lifetime(15) = 1.0E+30 ! Antiproton
+ geant_lifetime(16) = 8.922E-11 ! Kaon 0 short
+ geant_lifetime(17) = 5.53085E-19 ! Eta
+ geant_lifetime(18) = 2.632E-10 ! Lambda
+ geant_lifetime(19) = 7.99E-11 ! Sigma+
+ geant_lifetime(20) = 7.40E-20 ! Sigma0
+ geant_lifetime(21) = 1.479E-10 ! Sigma-
+ geant_lifetime(22) = 2.90E-10 ! Xi 0
+ geant_lifetime(23) = 1.639E-10 ! Xi -
+ geant_lifetime(24) = 8.22E-11 ! Omega
+ geant_lifetime(25) = 889.1 ! Antineutron
+ geant_lifetime(26) = 2.632E-10 ! Antilambda
+ geant_lifetime(27) = 7.99E-11 ! Anti-Sigma -
+ geant_lifetime(28) = 7.40E-20 ! Anti-Sigma 0
+ geant_lifetime(29) = 1.479E-10 ! Anti-Sigma +
+ geant_lifetime(30) = 2.90E-10 ! AntiXi 0
+ geant_lifetime(31) = 1.639E-10 ! AntiXi +
+ geant_lifetime(32) = 8.22E-11 ! Anti-Omega +
+ geant_lifetime(33) = 0.303E-12 ! Tau +
+ geant_lifetime(34) = 0.303E-12 ! Tau -
+ geant_lifetime(35) = 10.62E-13 ! D+
+ geant_lifetime(36) = 10.62E-13 ! D-
+ geant_lifetime(37) = 4.21E-13 ! D0
+ geant_lifetime(38) = 4.21E-13 ! Anti D0
+ geant_lifetime(39) = 4.45E-13 ! F+, now called Ds+
+ geant_lifetime(40) = 4.45E-13 ! F-, now called Ds-
+ geant_lifetime(41) = 1.91E-13 ! Lambda C+
+ geant_lifetime(42) = 2.92E-25 ! W+
+ geant_lifetime(43) = 2.92E-25 ! W-
+ geant_lifetime(44) = 2.60E-25 ! Z0
+ geant_lifetime(45) = 1.0E+30 ! Deuteron
+ geant_lifetime(46) = 1.0E+30 ! Triton
+ geant_lifetime(47) = 1.0E+30 ! Alpha
+ geant_lifetime(48) = 1.0E+30 ! Geantino (Fake particle)
+ geant_lifetime(49) = 1.0E+30 ! He3
+ geant_lifetime(50) = 1.0E+30 ! Cerenkov (Fake particle)
+ geant_lifetime(151) = 3.72E-24 ! rho +
+ geant_lifetime(152) = 3.72E-24 ! rho -
+ geant_lifetime(153) = 3.72E-24 ! rho 0
+ geant_lifetime(154) = 7.84E-23 ! omega 0
+ geant_lifetime(155) = 3.16E-21 ! eta'
+ geant_lifetime(156) = 1.49E-22 ! phi
+ geant_lifetime(157) = 9.68E-21 ! J/Psi
+ geant_lifetime(158) = 9.27E-24 ! Delta -, Based on gamma = 71 MeV
+ geant_lifetime(159) = 9.27E-24 ! Delta 0, -same-
+ geant_lifetime(160) = 9.27E-24 ! Delta +, -same-
+ geant_lifetime(161) = 9.27E-24 ! Delta ++,-same-
+ geant_lifetime(162) = 1.322E-23 ! K* +
+ geant_lifetime(163) = 1.322E-23 ! K* -
+ geant_lifetime(164) = 1.303E-23 ! K* 0
+
+CCC Set Particle Widths in GeV:
+ do i = 1,200
+ if(geant_lifetime(i) .le. 0.0) then
+ geant_width(i) = 0.0
+ else if(geant_lifetime(i) .ge. 1.0E+30) then
+ geant_width(i) = 0.0
+ else
+ geant_width(i) = hbar/geant_lifetime(i)
+ end if
+ end do
+
+ Return
+ END
+
+ Real*4 Function Vn_pt_y(n,V1,V2,V3,V4,pt,y)
+ implicit none
+
+CCC Description: This function computes the pt,y dependent flow
+CCC parameters Vn(pt,y) for the requested Fourier
+CCC component, n = 1-6, pt (GeV/c) and y (rapidity).
+CCC
+CCC Input: n = Fourier component, 1,2...6
+CCC V1 = Constant coefficient in pt dependent term
+CCC V2 = Coefficient of pt(pt^2) dependence for n odd (even).
+CCC V3 = Coefficient of y dependence; constant for n=odd,
+CCC and inverse range squared for Gaussian for n=even.
+CCC V4 = Coefficient of y^3 dependence for n=odd;
+CCC pt dependence for n = even.
+CCC pt = Transverse momentum (GeV/c)
+CCC y = Rapidity relative to y(C.M.)
+CCC
+CCC Output: Vn_pt_y = Vn(pt,y) based on the model suggested by
+CCC Art Poskanzer (LBNL, Feb. 2000)
+CCC Vn_pt_y = (V1 + V4*pt + V2*pt*pt)*exp(-V3*y*y) ; n=even
+CCC Vn_pt_y = (V1 + V2*pt)*sign(y)*(V3 + V4*abs(y**3)); n=odd
+
+CCC Local Variable Type Declarations:
+
+ integer n
+ real*4 V1,V2,V3,V4,pt,y,signy
+
+ if(n .eq. (2*(n/2))) then
+ Vn_pt_y = (V1 + V4*pt + V2*pt*pt)*exp(-V3*y*y)
+ else
+ if(y.ge.0.0) then
+ signy = 1.0
+ else if(y.lt.0.0) then
+ signy = -1.0
+ end if
+ Vn_pt_y = (V1 + V2*pt)*signy*(V3 + V4*abs(y**3))
+ end if
+ Return
+ END
+
+ Subroutine Particle_mass(gpid,pid_index,npts)
+ implicit none
+
+CCC Description: This subroutine computes the mass distributions for
+C included resonances at 'npts' number of mesh points in mass from the
+C lower mass limit to an upper mass limit (listed below), integrates
+C this mass distribution, normalizes the integral to 1.0, and saves
+C the mass steps and integral function in the arrays in Common/Mass/.
+C The mass distribution integral is then randomly sampled in a later
+C step in order to get the specific resonance mass instances.
+C For non-resonant particles (i.e. either stable or those with negligible
+C widths) this subroutine returns without doing anything, leaving the
+C arrays in Common/Mass/ set to zero. This subroutine is called for
+C a specific PID index, corresponding to the input list of particle
+C types.
+C
+C Input: gpid = Geant particle ID code number, see SUBR:
+C Particle_prop for listing.
+C pid_index = Particle type array index, determined by input
+C particle list.
+C npts = n_integ_pts in calling program; is the number
+C of mass mesh points used to load the mass
+C distribution integral. Note that one extra
+C mesh point is added to handle the bug in the
+C Lagrange interpolator, LAGRNG.
+C
+C Output: Mass_integ_save( , ) - mass distribution integral
+C Mass_xfunc_save( , ) - mass distribution mesh
+C These are in Common/Mass/.
+
+CCC Include files and common blocks:
+ Include 'Parameter_values.inc'
+ Common/Geant/geant_mass(200),geant_charge(200),
+ 1 geant_lifetime(200),geant_width(200)
+ real*4 geant_mass,geant_charge,geant_lifetime,geant_width
+ Common/Mass/ Mass_integ_save(npid,nmax_integ),
+ 1 Mass_xfunc_save(npid,nmax_integ)
+ real*4 Mass_integ_save,Mass_xfunc_save
+
+CCC Local Variable Type Declarations:
+ integer gpid,pid_index,npts,i
+ real*4 dist(nmax_integ),dM,M0,Gamma,GammaS
+ real*4 M,Mpi,MK,MN,R_Delta,Jinc,qR
+ real*4 Mrho_low,Mrho_high,Momega_low,Momega_high
+ real*4 Mphi_low,Mphi_high,MDelta_low,MDelta_high
+ real*4 MKstar_low,MKstar_high
+ real*4 kcm,k0cm,Ecm,E0cm,beta,beta0,Epicm,ENcm,redtotE
+
+CCC Set Fixed parameter values:
+ Parameter(Mpi = 0.1395675) ! Charged pion mass (GeV)
+ Parameter(MK = 0.493646) ! Charged kaon mass (GeV)
+ Parameter(MN = 0.938919) ! Nucleon average mass (GeV)
+ Parameter(R_Delta = 0.81) ! Delta resonance range parameter
+ Parameter(Mrho_low = 0.28 ) ! Lower rho mass limit
+ Parameter(Mrho_high = 1.200) ! Upper rho mass limit (GeV)
+ Parameter(Momega_low = 0.75) ! Lower omega mass limit (GeV)
+ Parameter(Momega_high = 0.81) ! Upper omega mass limit (GeV)
+ Parameter(Mphi_low = 1.009) ! Lower phi mass limit (GeV)
+ Parameter(Mphi_high = 1.029) ! Upper phi mass limit (GeV)
+ Parameter(MDelta_low = 1.1 ) ! Lower Delta mass limit
+ Parameter(MDelta_high = 1.400) ! Upper Delta mass limit (GeV)
+ Parameter(MKstar_low = 0.74) ! Lower Kstar mass limit (GeV)
+ Parameter(MKstar_high = 1.04) ! Upper Kstar mass limit (GeV)
+
+CCC Check npts:
+ if(npts.le.1) Return
+
+CCC Load mass distribution for rho-meson:
+ if(gpid.ge.151 .and. gpid.le.153) then
+ do i = 1,nmax_integ
+ dist(i) = 0.0
+ end do
+ M0 = geant_mass(gpid)
+ Gamma = geant_width(gpid)
+ dM = (Mrho_high - Mrho_low)/float(npts-1)
+ do i = 1,npts
+ M = Mrho_low + dM*float(i-1)
+ Mass_xfunc_save(pid_index,i) = M
+ kcm = sqrt(0.25*M*M - Mpi*Mpi)
+ dist(i) = kcm/((M-M0)**2 + 0.25*Gamma*Gamma)
+ end do
+
+CCC Load mass distribution for omega-meson:
+ else if(gpid.eq.154) then
+ do i = 1,nmax_integ
+ dist(i) = 0.0
+ end do
+ M0 = geant_mass(gpid)
+ Gamma = geant_width(gpid)
+ dM = (Momega_high - Momega_low)/float(npts-1)
+ do i = 1,npts
+ M = Momega_low + dM*float(i-1)
+ Mass_xfunc_save(pid_index,i) = M
+ GammaS = Gamma*((M/M0)**3)
+ dist(i) = M0*M0*Gamma/((M0*M0 - M*M)**2
+ 1 + M*M*GammaS*GammaS)
+ end do
+
+CCC Load mass distribution for phi-meson:
+ else if(gpid.eq.156) then
+ do i = 1,nmax_integ
+ dist(i) = 0.0
+ end do
+ M0 = geant_mass(gpid)
+ Gamma = geant_width(gpid)
+ dM = (Mphi_high - Mphi_low)/float(npts-1)
+ k0cm = sqrt(0.25*M0*M0 - MK*MK)
+ E0cm = sqrt(k0cm*k0cm + MK*MK)
+ beta0 = k0cm/E0cm
+ do i = 1,npts
+ M = Mphi_low + dM*float(i-1)
+ Mass_xfunc_save(pid_index,i) = M
+ kcm = sqrt(0.25*M*M - MK*MK)
+ Ecm = sqrt(kcm*kcm + MK*MK)
+ beta = kcm/Ecm
+ dist(i) = 0.25*Gamma*Gamma*((beta/beta0)**3)/
+ 1 ((M - M0)**2 + 0.25*Gamma*Gamma)
+ end do
+
+CCC Load mass distribution for Delta resonances:
+ else if(gpid.ge.158 .and. gpid.le.161) then
+ do i = 1,nmax_integ
+ dist(i) = 0.0
+ end do
+ M0 = geant_mass(gpid)
+ Gamma = geant_width(gpid)
+ dM = (MDelta_high - MDelta_low)/float(npts-1)
+ do i = 1,npts
+ M = MDelta_low + dM*float(i-1)
+ Mass_xfunc_save(pid_index,i) = M
+ kcm = ((M*M + Mpi*Mpi - MN*MN)**2)/(4.0*M*M) - Mpi*Mpi
+ kcm = sqrt(abs(kcm))
+ Epicm = sqrt(kcm*kcm + Mpi*Mpi)
+ ENcm = sqrt(kcm*kcm + MN*MN)
+ redtotE = Epicm*ENcm/(Epicm + ENcm)
+ Jinc = kcm/redtotE
+ qR = kcm*R_Delta/Mpi
+ GammaS = 2.0*qR*qR*qR*Gamma/(1.0 + qR*qR)
+ dist(i) = (Jinc/(kcm*kcm))*GammaS*GammaS/
+ 1 ((M - M0)**2 + 0.25*GammaS*GammaS)
+ end do
+
+CCC Load mass distribution for K*(892) resonances:
+ else if(gpid.ge.162 .and. gpid.le.164) then
+ do i = 1,nmax_integ
+ dist(i) = 0.0
+ end do
+ M0 = geant_mass(gpid)
+ Gamma = geant_width(gpid)
+ dM = (MKstar_high - MKstar_low)/float(npts-1)
+ k0cm = ((M0*M0 + Mpi*Mpi - MK*MK)**2)/(4.0*M0*M0) - Mpi*Mpi
+ k0cm = sqrt(k0cm)
+ do i = 1,npts
+ M = MKstar_low + dM*float(i-1)
+ Mass_xfunc_save(pid_index,i) = M
+ kcm = ((M*M + Mpi*Mpi - MK*MK)**2)/(4.0*M*M) - Mpi*Mpi
+ kcm = sqrt(kcm)
+ qR = kcm/k0cm
+ GammaS = 2.0*qR*qR*qR*Gamma/(1.0 + qR*qR)
+ dist(i) = GammaS*GammaS*M0*M0/
+ 1 ((M*M - M0*M0)**2 + GammaS*GammaS*M0*M0)
+ end do
+
+CCC Load additional resonance mass distributions here:
+
+ else
+ Return ! Return for Geant PID types without mass distribution
+ end if
+
+CCC Integrate mass distribution from mass_low to mass_high:
+
+ Mass_integ_save(pid_index,1) = 0.0
+ do i = 2,npts
+ Mass_integ_save(pid_index,i) = 0.5*(dist(i-1) + dist(i))*dM
+ 1 + Mass_integ_save(pid_index,i-1)
+ end do
+
+CCC Normalize this integral such that the last point is 1.00:
+ if(Mass_integ_save(pid_index,npts) .ne. 0.0) then
+ do i = 1,npts
+ Mass_integ_save(pid_index,i) =
+ 1 Mass_integ_save(pid_index,i)/Mass_integ_save(pid_index,npts)
+ end do
+ end if
+
+CCC Extend integ() and xfunc() by one mesh point to avoid Lagrange
+CCC interpolation subroutine bug:
+ Mass_integ_save(pid_index,npts+1) =
+ 1 Mass_integ_save(pid_index,npts) + 0.01
+ Mass_xfunc_save(pid_index,npts+1) =
+ 1 Mass_xfunc_save(pid_index,npts)
+
+ Return
+ END
+
+ Real*4 Function Mass_function(gpid,pid_index,npts,irand)
+ implicit none
+
+CCC Description: For resonance particles which have mass distributions
+C this function randomly samples the distributions in Common/Mass/
+C and returns a particle mass in GeV in 'Mass_function'.
+C For non-resonant particles this function returns the Geant mass
+C listed in SUBR: Particle_prop.
+C
+C Input: gpid = Geant particle ID code number, see SUBR:
+C Particle_prop for listing.
+C pid_index = Particle type array index, determined by input
+C particle list.
+C npts = n_integ_pts in calling program. Is the number
+C of mass mesh points for the arrays
+C in Common/Mass/ minus 1.
+C irand = random number generator argument/seed
+C
+C Output: Mass_function = particle or resonance mass (GeV)
+
+CCC Include files and common blocks:
+ Include 'Parameter_values.inc'
+ Common/Geant/geant_mass(200),geant_charge(200),
+ 1 geant_lifetime(200),geant_width(200)
+ real*4 geant_mass,geant_charge,geant_lifetime,geant_width
+ Common/Mass/ Mass_integ_save(npid,nmax_integ),
+ 1 Mass_xfunc_save(npid,nmax_integ)
+ real*4 Mass_integ_save,Mass_xfunc_save
+
+CCC Local Variable Type Declarations:
+ integer gpid,pid_index,npts,irand,i
+ real*4 integ(nmax_integ),xfunc(nmax_integ)
+ real*4 ran,masstmp
+
+ if(Mass_integ_save(pid_index,npts) .ne. 0.0) then
+ do i = 1,npts+1
+ integ(i) = Mass_integ_save(pid_index,i)
+ xfunc(i) = Mass_xfunc_save(pid_index,i)
+ end do
+ Call LAGRNG(ran(irand),integ,masstmp,xfunc,
+ 1 npts+1, 1, 5, npts+1, 1)
+ Mass_function = masstmp
+ else
+ Mass_function = geant_mass(gpid)
+ end if
+
+ Return
+ END
+*
+* real*4 function ran(i)
+* implicit none
+* integer i
+* real*4 r
+* Call ranlux2(r,1,i)
+* ran = r
+* Return
+* END
+*
+* Include 'ranlux2.F'
+
+
+
+
+
+
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff