From c28b8f8ddd03472896cd635bfc69a1d5c6f7e357 Mon Sep 17 00:00:00 2001 From: morsch Date: Sun, 25 Mar 2001 10:17:10 +0000 Subject: [PATCH] First commit of the MeVSim Fortran code. --- MEVSIM/Makefile | 70 + MEVSIM/Parameter_values.inc | 14 + MEVSIM/common_dist_bin.inc | 19 + MEVSIM/common_facfac.inc | 9 + MEVSIM/dummymevsim.F | 2 + MEVSIM/multiplicity_gen.F | 2755 +++++++++++++++++++++++++++++++++++ MEVSIM/ranlux.F | 309 ++++ MEVSIM/ranlux2.F | 316 ++++ 8 files changed, 3494 insertions(+) create mode 100644 MEVSIM/Makefile create mode 100644 MEVSIM/Parameter_values.inc create mode 100644 MEVSIM/common_dist_bin.inc create mode 100644 MEVSIM/common_facfac.inc create mode 100644 MEVSIM/dummymevsim.F create mode 100644 MEVSIM/multiplicity_gen.F create mode 100644 MEVSIM/ranlux.F create mode 100644 MEVSIM/ranlux2.F diff --git a/MEVSIM/Makefile b/MEVSIM/Makefile new file mode 100644 index 00000000000..0d9b0898620 --- /dev/null +++ b/MEVSIM/Makefile @@ -0,0 +1,70 @@ +############################### MEVSIM Makefile ############################### + +# Include machine specific definitions + +include $(ALICE_ROOT)/conf/GeneralDef +include $(ALICE_ROOT)/conf/MachineDef.$(ALICE_TARGET) + +PACKAGE = MEVSIM + + +# C++ sources + + +##### MACROS ##### + +FSRCS = multiplicity_gen.F + +FOBJS = $(patsubst %.F,tgt_$(ALICE_TARGET)/%.o,$(FSRCS)) + +SRCS = $(FSRCS) $(CSRCS) +OBJS = $(FOBJS) + +DSRCS = dummymevsim.F +DOBJS = $(patsubst %.F,tgt_$(ALICE_TARGET)/%.o,$(DSRCS)) + +# C++ compilation flags + +CXXFLAGS = $(CXXOPTS) $(CLIBCXXOPTS) $(CLIBDEFS) + + +# C compilation flags + +CFLAGS = $(COPT) $(CLIBCOPT) $(CLIBDEFS) + +# FORTRAN compilation flags + +FFLAGS = $(FOPT) $(CLIBFOPT) $(CLIBDEFS) + +##### TARGETS ##### + +# Target + +SLIBRARY = $(LIBDIR)/libmevsim.$(SL) $(LIBDIR)/libdummymevsim.$(SL) +ALIBRARY = $(LIBDIR)/libmevsim.a + +default: $(SLIBRARY) + +$(LIBDIR)/libmevsim.$(SL): $(OBJS) +$(LIBDIR)/libdummymevsim.$(SL): $(DOBJS) + +depend: $(CSRCS) $(DSRCS) $(SRCS) + +TOCLEAN = $(OBJS) $(DOBJS) *Cint.cxx *Cint.h + +############################### General Macros ################################ + +include $(ALICE_ROOT)/conf/GeneralMacros + +############################ Dependencies ##################################### + +-include tgt_$(ALICE_TARGET)/Make-depend + + + + + + + + + diff --git a/MEVSIM/Parameter_values.inc b/MEVSIM/Parameter_values.inc new file mode 100644 index 00000000000..3d4633e1a69 --- /dev/null +++ b/MEVSIM/Parameter_values.inc @@ -0,0 +1,14 @@ +CCC Set array dimension sizes here: + + integer npid,nmax_integ,n_mult_max_steps,nflowterms + parameter (npid = 30) ! max # of particle ID types + parameter (nmax_integ = 100) ! max # integration steps in parameter +CCC ! variance calculation. + parameter (n_mult_max_steps = 1000) +CCC ! max # integration steps in multiplicity +CCC ! variance calculation (this must be an +CCC ! even integer). + parameter (nflowterms = 6) ! max # of terms in the anisotropic +CCC ! flow model for azimuthal (phi angle) +CCC ! dependence. + diff --git a/MEVSIM/common_dist_bin.inc b/MEVSIM/common_dist_bin.inc new file mode 100644 index 00000000000..5ae91650443 --- /dev/null +++ b/MEVSIM/common_dist_bin.inc @@ -0,0 +1,19 @@ +CCC Common for bin-by-bin distribution input: +CCC NOTE: Include file 'Parameter_values.inc' must accompany and +CCC precede this file everywhere it occurs. + + integer n_bins_max + parameter(n_bins_max = 50) ! maximum # of input pt, eta bins + + Common/dist_bin/ pt_start(npid),eta_start(npid),pt_stop(npid), + 1 eta_stop(npid),delta_pt(npid,n_bins_max), + 2 delta_eta(npid,n_bins_max),pt_bin(npid,n_bins_max), + 3 eta_bin(npid,n_bins_max), + 4 pt_eta_bin(npid,n_bins_max,n_bins_max), + 5 pt_bin_mesh(npid,n_bins_max),eta_bin_mesh(npid,n_bins_max), + 6 n_pt_bins(npid),n_eta_bins(npid) + + integer n_pt_bins,n_eta_bins + real*4 pt_start,eta_start,pt_stop,eta_stop + real*4 delta_pt,delta_eta,pt_bin,eta_bin,pt_eta_bin + real*4 pt_bin_mesh,eta_bin_mesh diff --git a/MEVSIM/common_facfac.inc b/MEVSIM/common_facfac.inc new file mode 100644 index 00000000000..870b350ac43 --- /dev/null +++ b/MEVSIM/common_facfac.inc @@ -0,0 +1,9 @@ +CCC Common FACFAC for Factorials + + integer factorial_max + parameter (factorial_max = 10000) ! max # multiplicity per event; +CCC ! for any specific particle ID; +CCC ! also used for log(n!). + Common/FACFAC/ FACLOG(factorial_max) + real*4 FACLOG + diff --git a/MEVSIM/dummymevsim.F b/MEVSIM/dummymevsim.F new file mode 100644 index 00000000000..b16dfab88b8 --- /dev/null +++ b/MEVSIM/dummymevsim.F @@ -0,0 +1,2 @@ + subroutine track() + end diff --git a/MEVSIM/multiplicity_gen.F b/MEVSIM/multiplicity_gen.F new file mode 100644 index 00000000000..4c7e10bcb6d --- /dev/null +++ b/MEVSIM/multiplicity_gen.F @@ -0,0 +1,2755 @@ +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' + + + + + + diff --git a/MEVSIM/ranlux.F b/MEVSIM/ranlux.F new file mode 100644 index 00000000000..030639a1408 --- /dev/null +++ b/MEVSIM/ranlux.F @@ -0,0 +1,309 @@ +C* +C* $Id$ +C* +C* $Log$ +C* Revision 1.2 1997/09/22 13:45:47 mclareni +C* Correct error in initializing RANLUX by using RLUXIN with the output of +C* RLUXUT from a previous run. +C* +C* Revision 1.1.1.1 1996/04/01 15:02:55 mclareni +C* Mathlib gen +C* +C* +C#include "gen/pilot.h" + SUBROUTINE RANLUX(RVEC,LENV) +C Subtract-and-borrow random number generator proposed by +C Marsaglia and Zaman, implemented by F. James with the name +C RCARRY in 1991, and later improved by Martin Luescher +C in 1993 to produce "Luxury Pseudorandom Numbers". +C Fortran 77 coded by F. James, 1993 +C +C LUXURY LEVELS. +C ------ ------ The available luxury levels are: +C +C level 0 (p=24): equivalent to the original RCARRY of Marsaglia +C and Zaman, very long period, but fails many tests. +C level 1 (p=48): considerable improvement in quality over level 0, +C now passes the gap test, but still fails spectral test. +C level 2 (p=97): passes all known tests, but theoretically still +C defective. +C level 3 (p=223): DEFAULT VALUE. Any theoretically possible +C correlations have very small chance of being observed. +C level 4 (p=389): highest possible luxury, all 24 bits chaotic. +C +C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +C!!! Calling sequences for RANLUX: ++ +C!!! CALL RANLUX (RVEC, LEN) returns a vector RVEC of LEN ++ +C!!! 32-bit random floating point numbers between ++ +C!!! zero (not included) and one (also not incl.). ++ +C!!! CALL RLUXGO(LUX,INT,K1,K2) initializes the generator from ++ +C!!! one 32-bit integer INT and sets Luxury Level LUX ++ +C!!! which is integer between zero and MAXLEV, or if ++ +C!!! LUX .GT. 24, it sets p=LUX directly. K1 and K2 ++ +C!!! should be set to zero unless restarting at a break++ +C!!! point given by output of RLUXAT (see RLUXAT). ++ +C!!! CALL RLUXAT(LUX,INT,K1,K2) gets the values of four integers++ +C!!! which can be used to restart the RANLUX generator ++ +C!!! at the current point by calling RLUXGO. K1 and K2++ +C!!! specify how many numbers were generated since the ++ +C!!! initialization with LUX and INT. The restarting ++ +C!!! skips over K1+K2*E9 numbers, so it can be long.++ +C!!! A more efficient but less convenient way of restarting is by: ++ +C!!! CALL RLUXIN(ISVEC) restarts the generator from vector ++ +C!!! ISVEC of 25 32-bit integers (see RLUXUT) ++ +C!!! CALL RLUXUT(ISVEC) outputs the current values of the 25 ++ +C!!! 32-bit integer seeds, to be used for restarting ++ +C!!! ISVEC must be dimensioned 25 in the calling program ++ +C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + DIMENSION RVEC(LENV) + DIMENSION SEEDS(24), ISEEDS(24), ISDEXT(25) + PARAMETER (MAXLEV=4, LXDFLT=3) + DIMENSION NDSKIP(0:MAXLEV) + DIMENSION NEXT(24) + PARAMETER (TWOP12=4096., IGIGA=1000000000,JSDFLT=314159265) + PARAMETER (ITWO24=2**24, ICONS=2147483563) + SAVE NOTYET, I24, J24, CARRY, SEEDS, TWOM24, TWOM12, LUXLEV + SAVE NSKIP, NDSKIP, IN24, NEXT, KOUNT, MKOUNT, INSEED + INTEGER LUXLEV + LOGICAL NOTYET + DATA NOTYET, LUXLEV, IN24, KOUNT, MKOUNT /.TRUE., LXDFLT, 0,0,0/ + DATA I24,J24,CARRY/24,10,0./ +C default +C Luxury Level 0 1 2 *3* 4 + DATA NDSKIP/0, 24, 73, 199, 365 / +Corresponds to p=24 48 97 223 389 +C time factor 1 2 3 6 10 on slow workstation +C 1 1.5 2 3 5 on fast mainframe +C +C NOTYET is .TRUE. if no initialization has been performed yet. +C Default Initialization by Multiplicative Congruential + IF (NOTYET) THEN + NOTYET = .FALSE. + JSEED = JSDFLT + INSEED = JSEED + WRITE(6,'(A,I12)') ' RANLUX DEFAULT INITIALIZATION: ',JSEED + LUXLEV = LXDFLT + NSKIP = NDSKIP(LUXLEV) + LP = NSKIP + 24 + IN24 = 0 + KOUNT = 0 + MKOUNT = 0 + WRITE(6,'(A,I2,A,I4)') ' RANLUX DEFAULT LUXURY LEVEL = ', + + LUXLEV,' p =',LP + TWOM24 = 1. + DO 25 I= 1, 24 + TWOM24 = TWOM24 * 0.5 + K = JSEED/53668 + JSEED = 40014*(JSEED-K*53668) -K*12211 + IF (JSEED .LT. 0) JSEED = JSEED+ICONS + ISEEDS(I) = MOD(JSEED,ITWO24) + 25 CONTINUE + TWOM12 = TWOM24 * 4096. + DO 50 I= 1,24 + SEEDS(I) = REAL(ISEEDS(I))*TWOM24 + NEXT(I) = I-1 + 50 CONTINUE + NEXT(1) = 24 + I24 = 24 + J24 = 10 + CARRY = 0. + IF (SEEDS(24) .EQ. 0.) CARRY = TWOM24 + ENDIF +C +C The Generator proper: "Subtract-with-borrow", +C as proposed by Marsaglia and Zaman, +C Florida State University, March, 1989 +C + DO 100 IVEC= 1, LENV + UNI = SEEDS(J24) - SEEDS(I24) - CARRY + IF (UNI .LT. 0.) THEN + UNI = UNI + 1.0 + CARRY = TWOM24 + ELSE + CARRY = 0. + ENDIF + SEEDS(I24) = UNI + I24 = NEXT(I24) + J24 = NEXT(J24) + RVEC(IVEC) = UNI +C small numbers (with less than 12 "significant" bits) are "padded". + IF (UNI .LT. TWOM12) THEN + RVEC(IVEC) = RVEC(IVEC) + TWOM24*SEEDS(J24) +C and zero is forbidden in case someone takes a logarithm + IF (RVEC(IVEC) .EQ. 0.) RVEC(IVEC) = TWOM24*TWOM24 + ENDIF +C Skipping to luxury. As proposed by Martin Luscher. + IN24 = IN24 + 1 + IF (IN24 .EQ. 24) THEN + IN24 = 0 + KOUNT = KOUNT + NSKIP + DO 90 ISK= 1, NSKIP + UNI = SEEDS(J24) - SEEDS(I24) - CARRY + IF (UNI .LT. 0.) THEN + UNI = UNI + 1.0 + CARRY = TWOM24 + ELSE + CARRY = 0. + ENDIF + SEEDS(I24) = UNI + I24 = NEXT(I24) + J24 = NEXT(J24) + 90 CONTINUE + ENDIF + 100 CONTINUE + KOUNT = KOUNT + LENV + IF (KOUNT .GE. IGIGA) THEN + MKOUNT = MKOUNT + 1 + KOUNT = KOUNT - IGIGA + ENDIF + RETURN +C +C Entry to input and float integer seeds from previous run + ENTRY RLUXIN(ISDEXT) + NOTYET = .FALSE. + TWOM24 = 1. + DO 195 I= 1, 24 + NEXT(I) = I-1 + 195 TWOM24 = TWOM24 * 0.5 + NEXT(1) = 24 + TWOM12 = TWOM24 * 4096. + WRITE(6,'(A)') ' FULL INITIALIZATION OF RANLUX WITH 25 INTEGERS:' + WRITE(6,'(5X,5I12)') ISDEXT + DO 200 I= 1, 24 + SEEDS(I) = REAL(ISDEXT(I))*TWOM24 + 200 CONTINUE + CARRY = 0. + IF (ISDEXT(25) .LT. 0) CARRY = TWOM24 + ISD = IABS(ISDEXT(25)) + I24 = MOD(ISD,100) + ISD = ISD/100 + J24 = MOD(ISD,100) + ISD = ISD/100 + IN24 = MOD(ISD,100) + ISD = ISD/100 + LUXLEV = ISD + IF (LUXLEV .LE. MAXLEV) THEN + NSKIP = NDSKIP(LUXLEV) + WRITE (6,'(A,I2)') ' RANLUX LUXURY LEVEL SET BY RLUXIN TO: ', + + LUXLEV + ELSE IF (LUXLEV .GE. 24) THEN + NSKIP = LUXLEV - 24 + WRITE (6,'(A,I5)') ' RANLUX P-VALUE SET BY RLUXIN TO:',LUXLEV + ELSE + NSKIP = NDSKIP(MAXLEV) + WRITE (6,'(A,I5)') ' RANLUX ILLEGAL LUXURY RLUXIN: ',LUXLEV + LUXLEV = MAXLEV + ENDIF + INSEED = -1 + RETURN +C +C Entry to ouput seeds as integers + ENTRY RLUXUT(ISDEXT) + DO 300 I= 1, 24 + ISDEXT(I) = INT(SEEDS(I)*TWOP12*TWOP12) + 300 CONTINUE + ISDEXT(25) = I24 + 100*J24 + 10000*IN24 + 1000000*LUXLEV + IF (CARRY .GT. 0.) ISDEXT(25) = -ISDEXT(25) + RETURN +C +C Entry to output the "convenient" restart point + ENTRY RLUXAT(LOUT,INOUT,K1,K2) + LOUT = LUXLEV + INOUT = INSEED + K1 = KOUNT + K2 = MKOUNT + RETURN +C +C Entry to initialize from one or three integers + ENTRY RLUXGO(LUX,INS,K1,K2) + IF (LUX .LT. 0) THEN + LUXLEV = LXDFLT + ELSE IF (LUX .LE. MAXLEV) THEN + LUXLEV = LUX + ELSE IF (LUX .LT. 24 .OR. LUX .GT. 2000) THEN + LUXLEV = MAXLEV + WRITE (6,'(A,I7)') ' RANLUX ILLEGAL LUXURY RLUXGO: ',LUX + ELSE + LUXLEV = LUX + DO 310 ILX= 0, MAXLEV + IF (LUX .EQ. NDSKIP(ILX)+24) LUXLEV = ILX + 310 CONTINUE + ENDIF + IF (LUXLEV .LE. MAXLEV) THEN + NSKIP = NDSKIP(LUXLEV) + WRITE(6,'(A,I2,A,I4)') ' RANLUX LUXURY LEVEL SET BY RLUXGO :', + + LUXLEV,' P=', NSKIP+24 + ELSE + NSKIP = LUXLEV - 24 + WRITE (6,'(A,I5)') ' RANLUX P-VALUE SET BY RLUXGO TO:',LUXLEV + ENDIF + IN24 = 0 + IF (INS .LT. 0) WRITE (6,'(A)') + + ' Illegal initialization by RLUXGO, negative input seed' + IF (INS .GT. 0) THEN + JSEED = INS + WRITE(6,'(A,3I12)') ' RANLUX INITIALIZED BY RLUXGO FROM SEEDS', + + JSEED, K1,K2 + ELSE + JSEED = JSDFLT + WRITE(6,'(A)')' RANLUX INITIALIZED BY RLUXGO FROM DEFAULT SEED' + ENDIF + INSEED = JSEED + NOTYET = .FALSE. + TWOM24 = 1. + DO 325 I= 1, 24 + TWOM24 = TWOM24 * 0.5 + K = JSEED/53668 + JSEED = 40014*(JSEED-K*53668) -K*12211 + IF (JSEED .LT. 0) JSEED = JSEED+ICONS + ISEEDS(I) = MOD(JSEED,ITWO24) + 325 CONTINUE + TWOM12 = TWOM24 * 4096. + DO 350 I= 1,24 + SEEDS(I) = REAL(ISEEDS(I))*TWOM24 + NEXT(I) = I-1 + 350 CONTINUE + NEXT(1) = 24 + I24 = 24 + J24 = 10 + CARRY = 0. + IF (SEEDS(24) .EQ. 0.) CARRY = TWOM24 +C If restarting at a break point, skip K1 + IGIGA*K2 +C Note that this is the number of numbers delivered to +C the user PLUS the number skipped (if luxury .GT. 0). + KOUNT = K1 + MKOUNT = K2 + IF (K1+K2 .NE. 0) THEN + DO 500 IOUTER= 1, K2+1 + INNER = IGIGA + IF (IOUTER .EQ. K2+1) INNER = K1 + DO 450 ISK= 1, INNER + UNI = SEEDS(J24) - SEEDS(I24) - CARRY + IF (UNI .LT. 0.) THEN + UNI = UNI + 1.0 + CARRY = TWOM24 + ELSE + CARRY = 0. + ENDIF + SEEDS(I24) = UNI + I24 = NEXT(I24) + J24 = NEXT(J24) + 450 CONTINUE + 500 CONTINUE +C Get the right value of IN24 by direct calculation + IN24 = MOD(KOUNT, NSKIP+24) + IF (MKOUNT .GT. 0) THEN + IZIP = MOD(IGIGA, NSKIP+24) + IZIP2 = MKOUNT*IZIP + IN24 + IN24 = MOD(IZIP2, NSKIP+24) + ENDIF +C Now IN24 had better be between zero and 23 inclusive + IF (IN24 .GT. 23) THEN + WRITE (6,'(A/A,3I11,A,I5)') + + ' Error in RESTARTING with RLUXGO:',' The values', INS, + + K1, K2, ' cannot occur at luxury level', LUXLEV + IN24 = 0 + ENDIF + ENDIF + RETURN + END diff --git a/MEVSIM/ranlux2.F b/MEVSIM/ranlux2.F new file mode 100644 index 00000000000..65a8c847932 --- /dev/null +++ b/MEVSIM/ranlux2.F @@ -0,0 +1,316 @@ +C* +C* $Id$ +C* +C* $Log$ +C* Revision 1.2 1997/09/22 13:45:47 mclareni +C* Correct error in initializing RANLUX by using RLUXIN with the output of +C* RLUXUT from a previous run. +C* +C* Revision 1.1.1.1 1996/04/01 15:02:55 mclareni +C* Mathlib gen +C* +C* +C#include "gen/pilot.h" + SUBROUTINE RANLUX2(RVEC,LENV,Input_seed) +C Subtract-and-borrow random number generator proposed by +C Marsaglia and Zaman, implemented by F. James with the name +C RCARRY in 1991, and later improved by Martin Luescher +C in 1993 to produce "Luxury Pseudorandom Numbers". +C Fortran 77 coded by F. James, 1993 +C +C LUXURY LEVELS. +C ------ ------ The available luxury levels are: +C +C level 0 (p=24): equivalent to the original RCARRY of Marsaglia +C and Zaman, very long period, but fails many tests. +C level 1 (p=48): considerable improvement in quality over level 0, +C now passes the gap test, but still fails spectral test. +C level 2 (p=97): passes all known tests, but theoretically still +C defective. +C level 3 (p=223): DEFAULT VALUE. Any theoretically possible +C correlations have very small chance of being observed. +C level 4 (p=389): highest possible luxury, all 24 bits chaotic. +C +C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +C!!! Calling sequences for RANLUX: ++ +C!!! CALL RANLUX (RVEC, LEN) returns a vector RVEC of LEN ++ +C!!! 32-bit random floating point numbers between ++ +C!!! zero (not included) and one (also not incl.). ++ +C!!! CALL RLUXGO(LUX,INT,K1,K2) initializes the generator from ++ +C!!! one 32-bit integer INT and sets Luxury Level LUX ++ +C!!! which is integer between zero and MAXLEV, or if ++ +C!!! LUX .GT. 24, it sets p=LUX directly. K1 and K2 ++ +C!!! should be set to zero unless restarting at a break++ +C!!! point given by output of RLUXAT (see RLUXAT). ++ +C!!! CALL RLUXAT(LUX,INT,K1,K2) gets the values of four integers++ +C!!! which can be used to restart the RANLUX generator ++ +C!!! at the current point by calling RLUXGO. K1 and K2++ +C!!! specify how many numbers were generated since the ++ +C!!! initialization with LUX and INT. The restarting ++ +C!!! skips over K1+K2*E9 numbers, so it can be long.++ +C!!! A more efficient but less convenient way of restarting is by: ++ +C!!! CALL RLUXIN(ISVEC) restarts the generator from vector ++ +C!!! ISVEC of 25 32-bit integers (see RLUXUT) ++ +C!!! CALL RLUXUT(ISVEC) outputs the current values of the 25 ++ +C!!! 32-bit integer seeds, to be used for restarting ++ +C!!! ISVEC must be dimensioned 25 in the calling program ++ +C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + DIMENSION RVEC(LENV) + DIMENSION SEEDS(24), ISEEDS(24), ISDEXT(25) + PARAMETER (MAXLEV=4, LXDFLT=3) + DIMENSION NDSKIP(0:MAXLEV) + DIMENSION NEXT(24) + PARAMETER (TWOP12=4096., IGIGA=1000000000,JSDFLT=314159265) + PARAMETER (ITWO24=2**24, ICONS=2147483563) + SAVE NOTYET, I24, J24, CARRY, SEEDS, TWOM24, TWOM12, LUXLEV + SAVE NSKIP, NDSKIP, IN24, NEXT, KOUNT, MKOUNT, INSEED + INTEGER LUXLEV + Integer Input_seed,JSDFLT_set + LOGICAL NOTYET + DATA NOTYET, LUXLEV, IN24, KOUNT, MKOUNT /.TRUE., LXDFLT, 0,0,0/ + DATA I24,J24,CARRY/24,10,0./ +CCC Set starting seed value: + If(Input_seed.gt.0) then + JSDFLT_set = Input_seed + Else + JSDFLT_set = JSDFLT + End If +C default +C Luxury Level 0 1 2 *3* 4 + DATA NDSKIP/0, 24, 73, 199, 365 / +Corresponds to p=24 48 97 223 389 +C time factor 1 2 3 6 10 on slow workstation +C 1 1.5 2 3 5 on fast mainframe +C +C NOTYET is .TRUE. if no initialization has been performed yet. +C Default Initialization by Multiplicative Congruential + IF (NOTYET) THEN + NOTYET = .FALSE. + JSEED = JSDFLT_set + INSEED = JSEED + WRITE(6,'(A,I12)') ' RANLUX DEFAULT INITIALIZATION: ',JSEED + LUXLEV = LXDFLT + NSKIP = NDSKIP(LUXLEV) + LP = NSKIP + 24 + IN24 = 0 + KOUNT = 0 + MKOUNT = 0 + WRITE(6,'(A,I2,A,I4)') ' RANLUX DEFAULT LUXURY LEVEL = ', + + LUXLEV,' p =',LP + TWOM24 = 1. + DO 25 I= 1, 24 + TWOM24 = TWOM24 * 0.5 + K = JSEED/53668 + JSEED = 40014*(JSEED-K*53668) -K*12211 + IF (JSEED .LT. 0) JSEED = JSEED+ICONS + ISEEDS(I) = MOD(JSEED,ITWO24) + 25 CONTINUE + TWOM12 = TWOM24 * 4096. + DO 50 I= 1,24 + SEEDS(I) = REAL(ISEEDS(I))*TWOM24 + NEXT(I) = I-1 + 50 CONTINUE + NEXT(1) = 24 + I24 = 24 + J24 = 10 + CARRY = 0. + IF (SEEDS(24) .EQ. 0.) CARRY = TWOM24 + ENDIF +C +C The Generator proper: "Subtract-with-borrow", +C as proposed by Marsaglia and Zaman, +C Florida State University, March, 1989 +C + DO 100 IVEC= 1, LENV + UNI = SEEDS(J24) - SEEDS(I24) - CARRY + IF (UNI .LT. 0.) THEN + UNI = UNI + 1.0 + CARRY = TWOM24 + ELSE + CARRY = 0. + ENDIF + SEEDS(I24) = UNI + I24 = NEXT(I24) + J24 = NEXT(J24) + RVEC(IVEC) = UNI +C small numbers (with less than 12 "significant" bits) are "padded". + IF (UNI .LT. TWOM12) THEN + RVEC(IVEC) = RVEC(IVEC) + TWOM24*SEEDS(J24) +C and zero is forbidden in case someone takes a logarithm + IF (RVEC(IVEC) .EQ. 0.) RVEC(IVEC) = TWOM24*TWOM24 + ENDIF +C Skipping to luxury. As proposed by Martin Luscher. + IN24 = IN24 + 1 + IF (IN24 .EQ. 24) THEN + IN24 = 0 + KOUNT = KOUNT + NSKIP + DO 90 ISK= 1, NSKIP + UNI = SEEDS(J24) - SEEDS(I24) - CARRY + IF (UNI .LT. 0.) THEN + UNI = UNI + 1.0 + CARRY = TWOM24 + ELSE + CARRY = 0. + ENDIF + SEEDS(I24) = UNI + I24 = NEXT(I24) + J24 = NEXT(J24) + 90 CONTINUE + ENDIF + 100 CONTINUE + KOUNT = KOUNT + LENV + IF (KOUNT .GE. IGIGA) THEN + MKOUNT = MKOUNT + 1 + KOUNT = KOUNT - IGIGA + ENDIF + RETURN +C +C Entry to input and float integer seeds from previous run + ENTRY RLUXIN(ISDEXT) + NOTYET = .FALSE. + TWOM24 = 1. + DO 195 I= 1, 24 + NEXT(I) = I-1 + 195 TWOM24 = TWOM24 * 0.5 + NEXT(1) = 24 + TWOM12 = TWOM24 * 4096. + WRITE(6,'(A)') ' FULL INITIALIZATION OF RANLUX WITH 25 INTEGERS:' + WRITE(6,'(5X,5I12)') ISDEXT + DO 200 I= 1, 24 + SEEDS(I) = REAL(ISDEXT(I))*TWOM24 + 200 CONTINUE + CARRY = 0. + IF (ISDEXT(25) .LT. 0) CARRY = TWOM24 + ISD = IABS(ISDEXT(25)) + I24 = MOD(ISD,100) + ISD = ISD/100 + J24 = MOD(ISD,100) + ISD = ISD/100 + IN24 = MOD(ISD,100) + ISD = ISD/100 + LUXLEV = ISD + IF (LUXLEV .LE. MAXLEV) THEN + NSKIP = NDSKIP(LUXLEV) + WRITE (6,'(A,I2)') ' RANLUX LUXURY LEVEL SET BY RLUXIN TO: ', + + LUXLEV + ELSE IF (LUXLEV .GE. 24) THEN + NSKIP = LUXLEV - 24 + WRITE (6,'(A,I5)') ' RANLUX P-VALUE SET BY RLUXIN TO:',LUXLEV + ELSE + NSKIP = NDSKIP(MAXLEV) + WRITE (6,'(A,I5)') ' RANLUX ILLEGAL LUXURY RLUXIN: ',LUXLEV + LUXLEV = MAXLEV + ENDIF + INSEED = -1 + RETURN +C +C Entry to ouput seeds as integers + ENTRY RLUXUT(ISDEXT) + DO 300 I= 1, 24 + ISDEXT(I) = INT(SEEDS(I)*TWOP12*TWOP12) + 300 CONTINUE + ISDEXT(25) = I24 + 100*J24 + 10000*IN24 + 1000000*LUXLEV + IF (CARRY .GT. 0.) ISDEXT(25) = -ISDEXT(25) + RETURN +C +C Entry to output the "convenient" restart point + ENTRY RLUXAT(LOUT,INOUT,K1,K2) + LOUT = LUXLEV + INOUT = INSEED + K1 = KOUNT + K2 = MKOUNT + RETURN +C +C Entry to initialize from one or three integers + ENTRY RLUXGO(LUX,INS,K1,K2) + IF (LUX .LT. 0) THEN + LUXLEV = LXDFLT + ELSE IF (LUX .LE. MAXLEV) THEN + LUXLEV = LUX + ELSE IF (LUX .LT. 24 .OR. LUX .GT. 2000) THEN + LUXLEV = MAXLEV + WRITE (6,'(A,I7)') ' RANLUX ILLEGAL LUXURY RLUXGO: ',LUX + ELSE + LUXLEV = LUX + DO 310 ILX= 0, MAXLEV + IF (LUX .EQ. NDSKIP(ILX)+24) LUXLEV = ILX + 310 CONTINUE + ENDIF + IF (LUXLEV .LE. MAXLEV) THEN + NSKIP = NDSKIP(LUXLEV) + WRITE(6,'(A,I2,A,I4)') ' RANLUX LUXURY LEVEL SET BY RLUXGO :', + + LUXLEV,' P=', NSKIP+24 + ELSE + NSKIP = LUXLEV - 24 + WRITE (6,'(A,I5)') ' RANLUX P-VALUE SET BY RLUXGO TO:',LUXLEV + ENDIF + IN24 = 0 + IF (INS .LT. 0) WRITE (6,'(A)') + + ' Illegal initialization by RLUXGO, negative input seed' + IF (INS .GT. 0) THEN + JSEED = INS + WRITE(6,'(A,3I12)') ' RANLUX INITIALIZED BY RLUXGO FROM SEEDS', + + JSEED, K1,K2 + ELSE + JSEED = JSDFLT_set + WRITE(6,'(A)')' RANLUX INITIALIZED BY RLUXGO FROM DEFAULT SEED' + ENDIF + INSEED = JSEED + NOTYET = .FALSE. + TWOM24 = 1. + DO 325 I= 1, 24 + TWOM24 = TWOM24 * 0.5 + K = JSEED/53668 + JSEED = 40014*(JSEED-K*53668) -K*12211 + IF (JSEED .LT. 0) JSEED = JSEED+ICONS + ISEEDS(I) = MOD(JSEED,ITWO24) + 325 CONTINUE + TWOM12 = TWOM24 * 4096. + DO 350 I= 1,24 + SEEDS(I) = REAL(ISEEDS(I))*TWOM24 + NEXT(I) = I-1 + 350 CONTINUE + NEXT(1) = 24 + I24 = 24 + J24 = 10 + CARRY = 0. + IF (SEEDS(24) .EQ. 0.) CARRY = TWOM24 +C If restarting at a break point, skip K1 + IGIGA*K2 +C Note that this is the number of numbers delivered to +C the user PLUS the number skipped (if luxury .GT. 0). + KOUNT = K1 + MKOUNT = K2 + IF (K1+K2 .NE. 0) THEN + DO 500 IOUTER= 1, K2+1 + INNER = IGIGA + IF (IOUTER .EQ. K2+1) INNER = K1 + DO 450 ISK= 1, INNER + UNI = SEEDS(J24) - SEEDS(I24) - CARRY + IF (UNI .LT. 0.) THEN + UNI = UNI + 1.0 + CARRY = TWOM24 + ELSE + CARRY = 0. + ENDIF + SEEDS(I24) = UNI + I24 = NEXT(I24) + J24 = NEXT(J24) + 450 CONTINUE + 500 CONTINUE +C Get the right value of IN24 by direct calculation + IN24 = MOD(KOUNT, NSKIP+24) + IF (MKOUNT .GT. 0) THEN + IZIP = MOD(IGIGA, NSKIP+24) + IZIP2 = MKOUNT*IZIP + IN24 + IN24 = MOD(IZIP2, NSKIP+24) + ENDIF +C Now IN24 had better be between zero and 23 inclusive + IF (IN24 .GT. 23) THEN + WRITE (6,'(A/A,3I11,A,I5)') + + ' Error in RESTARTING with RLUXGO:',' The values', INS, + + K1, K2, ' cannot occur at luxury level', LUXLEV + IN24 = 0 + ENDIF + ENDIF + RETURN + END -- 2.39.3