Part of Reset() moved to base class.
[u/mrichter/AliRoot.git] / TMEVSIM / TMevSim.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 ////////////////////////////////////////////////////////////////////////////
19 // 
20 // TMevSim
21 //
22 // TMevSim is an interface class between the event generator MEVSIM and 
23 // the ROOT system. The current implementation is based on the 6.11.2000
24 // version provided by Lanny Ray on /afs/cern.ch/user/y/yiota/ray/mult_gen.
25 // 
26 // Authors of MEVSIM:
27 // For The STAR Collaboration
28 //
29 //  Lanny Ray      
30 //     Dept. of Physics
31 //     The University of Texas at Austin
32 //     Austin, Texas 78712
33 //     (512) 471-6107
34 //     ray@physics.utexas.edu
35 //
36 //  Ron Longacre   email:
37 //
38 //
39 ////////////////////////////////////////////////////////////////////////////
40 //
41 // I.    OVERVIEW
42 //
43 //    This code is intended to provide a quick means of producing
44 //    uncorrelated simulated events for event-by-event studies,
45 //    detector acceptance and efficiency studies, etc.  The
46 //    user selects the number of events, the one-particle distribution
47 //    model, the Geant particles to include, the ranges in transverse
48 //    momentum, pseudorapidity and azimuthal angle, the mean
49 //    multiplicity for each particle type for the event run, the
50 //    mean temperature, Rapidity width, etc., and the standard deviations
51 //    for the event-to-event variation in the model parameters.
52 //    Note that these events are produced in the c.m. frame only.
53 //    
54 //    Anisotropic flow may also be simulated by introducing explicit
55 //    phi-dependence (azimuthal angle) in the particle distributions.  
56 //    The assumed model is taken from Poskanzer and Voloshin, Phys. Rev.
57 //    C58, 1671 (1998), Eq.(1), where we use,
58 //
59 //         E d^3N/dp^3 = (1/2*pi*pt)*[d^2N/dpt*dy]
60 //            * [1 + SUM(n=1,nflowterms){2*Vn*cos[n(phi-PSIr)]}]
61 //
62 //    with up to 'nflowterms' (currently set to 6, see file
63 //    Parameter_values.inc) Fourier components allowed.  Vn are
64 //    coefficients and PSIr is the reaction plane angle.
65 //    The algebraic signs of the Vn terms for n=odd are reversed
66 //    from their input values for particles with rapidity (y) < 0
67 //    as suggested in Poskanzer and Voloshin.
68 //    The flow parameters can depend on pt and rapidity (y) according
69 //    to the model suggested by Art Poskanzer (Feb. 2000) and as
70 //    defined in the Function Vn_pt_y.
71 //
72 //    The user may select either to have the same multiplicity per
73 //    particle type for each event or to let the multiplicity vary
74 //    randomly according to a Poisson distribution. In addition, an
75 //    overall multiplicative scale factor can be applied to each
76 //    particle ID's multiplicity (same factor applied to each PID).
77 //    This scaling can vary randomly according to a Gaussian from
78 //    event-to-event.  This is to simulate trigger acceptance
79 //    fluctuations.  Similarly the
80 //    parameters of the one-particle distribution models may either
81 //    be fixed to the same value for each event or allowed to randomly
82 //    vary about a specified mean with a specified standard deviation
83 //    and assuming a Gaussian distribution.
84 //
85 //    With respect to the reaction plane and anisotropic flow simulation,
86 //    the user may select among four options:
87 //       (1) ignore reaction plane and anisotropic flow effects; all
88 //           distributions will be azimuthally invariant, on average.
89 //       (2) assume a fixed reaction plane angle, PSIr, for all events
90 //           in the run.
91 //       (3) assume a Gaussian distribution with specified mean and
92 //           standard deviation for the reaction plane angles for the
93 //           events in the run.  PSIr is randomly determined for each
94 //           event.
95 //       (4) assume uniformly distributed, random values for the reaction  
96 //           plane angles from 0 to 360 deg., for each event in the run.
97 //
98 //    The user may also select the anisotropic flow parameters, Vn,
99 //    to either be fixed for each event, or to randomly vary from event
100 //    to event according to a Gaussian distribution where the user must
101 //    specify the mean and std. dev.  For both cases the input file must
102 //    list the 'nflowterms' (e.g. 6) values of the mean and Std.dev. for
103 //    the Vn parameters for all particle ID types included in the run.
104 //
105 //    The available list of particles has been increased to permit a
106 //    number of meson and baryon resonances.  For those with broad widths
107 //    the code samples the mass distribution for the resonance and outputs
108 //    the resonance mass for each instance in a special kinematic file
109 //    (see file unit=9, filename = 'mult_gen.kin').  The resonance shapes
110 //    are approximately Breit-Wigner and are specific for each resonance 
111 //    case.  The additional particle/resonances include: rho(+,-,0),
112 //    omega(0), eta', phi, J/Psi, Delta(-,0,+,++) and K*(+,-,0).  Masses
113 //    are sampled for the rho, omega, phi, Deltas and D*s. 
114 //    Refer to SUBR: Particle_prop and Particle_mass for the explicit
115 //    parameters, resonance shape models, and sampling ranges.
116 //
117 //    The input is from a file, named 'mult_gen.in'.  The output is
118 //    loaded into a file named 'mult_gen.out' which includes run
119 //    header information, event header information and the EVENT: and
120 //    TRACK: formats as in the new STAR TEXT Format for event generator
121 //    input to GSTAR.  A log file, 'mult_gen.log' is also written which
122 //    may contain error messages.  Normally this file should be empty
123 //    after a successful run.  These filenames can easily be changed
124 //    to more suitable names by the script that runs the program or
125 //    by hand.
126 //
127 //
128 // II.   ALGORITHM
129 //
130 //
131 //
132 //    The method for generating random multiplicities and model parameter
133 //    values involves the following steps:
134 //       (1) The Poisson or Gaussian distributions are computed and
135 //           loaded into function f().
136 //       (2) The distribution f(x') is integrated from xmin to x
137 //           and saved from x = xmin to x = xmax.  The range and mesh
138 //           spaces are specified by the user.
139 //       (3) The integral of f is normalized to unity where 
140 //           integral[f(x')](at x = xmin) = 0.0
141 //           integral[f(x')](at x = xmax) = 1.0
142 //       (4) A random number generator is called which delivers values
143 //           between 0.0 and 1.0.  
144 //       (5) We consider the coordinate x (from xmin to xmax) to be
145 //           dependent on the integral[f].  Using the random number
146 //           for the selected value of integral[f] the value of x
147 //           is obtained by interpolation.
148 //
149 //    An interpolation subroutine from Rubin Landau, Oregon State Univ.,
150 //    is used to do this interpolation; it involves uneven mesh point 
151 //    spacing.
152 //
153 //    The method for generating the particle momenta uses the
154 //    standard random elimination method and involves the following
155 //    steps:
156 //
157 //    For model_type = 1,2,3,4 which are functions of pt,y (see following):
158 //       (1) The y range is computed using the pseudorapidity (eta)
159 //           range and includes ample cushioning around the sides
160 //           along the eta acceptance edges.
161 //       (2) The transverse momentum (pt) and rapidity (y) are
162 //           randomly chosen within the specified ranges.
163 //       (3) The pseudorapidity is computed for this (pt,y) value
164 //           (and the mass for each pid) and checked against the
165 //           pseudorapidity acceptance range.
166 //       (4) If the pseudorapidity is within range then the one-particle
167 //           model distribution is calculated at this point and its ratio
168 //           to the maximum value throughout (pt,eta) acceptance region
169 //           is calculated.
170 //       (5) Another random number is called and if less than the ratio
171 //           from step#4 the particle momentum is used; if not, then 
172 //           another trial value of (pt,y) is obtained.
173 //       (6) This continues until the required multiplicity for the
174 //           specific event and particle type has been satisfied.
175 //       (7) This process is repeated for the requested number of particle
176 //           types and events.
177 //
178 //    For model_type = 5,6 (see following) which are input bin-by-bin
179 //    in pt,eta:
180 //       (1) The transverse momentum (pt) and pseudorapidity (eta) are 
181 //           randomly chosen within the specified ranges.
182 //       (2) The one-particle model distribution is calculated at this
183 //           point and its ratio to the maximum value throughout the
184 //           (pt,eta) region is calculated.
185 //       (3) Another random number is called and if less than the ratio
186 //           from step(2) the particle momentum is used; if not then
187 //           another trial value of (pt,eta) is obtained.
188 //       (4) This continues until the required multiplicity for the 
189 //           specific event and particle type has been satisfied.
190 //       (5) This process is repeated for the requested number of particle
191 //           types and events. 
192 //
193 //    Problematic parameter values are tested, bad input values are checked
194 //    and in some cases may be changed so that the program will not crash.
195 //    In some cases the code execution is stopped.
196 //    Some distributions and/or unusual model parameter values may cause the
197 //    code to hang up due to the poor performance of the "elimination"
198 //    method for very strongly peaked distributions.  These are tested for
199 //    certain problematic values and if necessary these events are aborted.
200 //    A message, "*** Event No.    2903 ABORTED:" for example is printed
201 //    in the 'mult_gen.out' file.  Temperatures .le. 0.01 GeV and rapidity
202 //    width parameters .le. 0.01 will cause the event to abort.
203 //
204 //
205 //
206 // III.  DESCRIPTION OF THE INPUT:
207 //
208 //
209 //    The input is described below in the 'read' statements and also in
210 //    the sample input file.  Some additional comments are as follows:
211 //
212 //    (1) n_events - Selected number of events in run. Can be anything
213 //                   .ge. 1.
214 //    (2) n_pid_type - Number of particle ID types to include in the
215 //                     particle list. e.g. pi(+) and pi(-) are counted
216 //                     separately.  The limit is set by parameter npid
217 //                     in the accompanying include file 'Parameter_values.inc'
218 //                     and is presently set at 20.
219 //    (3) model_type - equals 1,2,3,4,5 or 6 so far.  See comments in
220 //                     Function dNdpty to see what is calculated.
221 //                     The models included are:
222 //                   = 1, Factorized mt exponential, Gaussian rapidity model
223 //                   = 2, Pratt non-expanding, spherical thermal source model
224 //                   = 3, Bertsch non-expanding spherical thermal source model
225 //                   = 4, Pratt spherically expanding, thermally equilibrated
226 //                        source model.
227 //                   = 5, Factorized pt and eta distributions input bin-by-bin.
228 //                   = 6, Fully 2D pt,eta distributions input bin-by-bin.
229 //                        NOTE: model_type = 1-4 are functions of (pt,y)
230 //                              model_type = 5,6 are functions of (pt,eta)
231 //    (4) reac_plane_cntrl - Can be either 1,2,3 or 4 where:
232 //                         = 1 to ignore reaction plane and anisotropic flow,
233 //                             all distributions will be azimuthally symm.
234 //                         = 2 to use a fixed reaction plane angle for all
235 //                             events in the run.
236 //                         = 3 to assume a randomly varying reaction plane
237 //                             angle for each event as determined by a
238 //                             Gaussian distribution.
239 //                         = 4 to assume a randomly varying reaction plane
240 //                             for each event in the run as determined by
241 //                             a uniform distribution from 0 to 360 deg.
242 //    (5) PSIr_mean, PSIr_stdev - Reaction plane angle mean and Gaussian
243 //                                std.dev. (both are in degrees) for cases
244 //                                with reac_plane_cntrl = 2 (use mean value)
245 //                                and 3.  Note: these are read in regardless
246 //                                of the value of reac_plane_cntrl.
247 //    (6) MultFac_mean, MultFac_stdev - Overall multiplicity scaling factor
248 //                                      for all PID types; mean and std.dev.;
249 //                                      for trigger fluctuations event-to-evt.
250 //    (7) pt_cut_min,pt_cut_max - Range of transverse momentum in GeV/c.
251 //    (8) eta_cut_min,eta_cut_max - Pseudorapidity range
252 //    (9) phi_cut_min,phi_cut_max - Azimuthal angular range in degrees.
253 //   (10) n_stdev_mult - Number of standard deviations about the mean value
254 //                       of multiplicity to include in the random event-to-
255 //                       event selection process.  The maximum number of
256 //                       steps that can be covered is determined by
257 //                       parameter n_mult_max_steps in the accompanying
258 //                       include file 'Parameter_values.inc' which is
259 //                       presently set at 1000, but the true upper limit for
260 //                       this is n_mult_max_steps - 1 = 999.
261 //   (11) n_stdev_temp - Same, except for the "Temperature" parameter.
262 //   (12) n_stdev_sigma- Same, except for the rapidity width parameter.
263 //   (13) n_stdev_expvel - Same, except for the expansion velocity parameter.
264 //   (14) n_stdev_PSIr   - Same, except for the reaction plane angle
265 //   (15) n_stdev_Vn     - Same, except for the anisotropy coefficients, Vn.
266 //   (16) n_stdev_MultFac - Same, except for the multiplicity scaling factor.
267 //   (17) n_integ_pts - Number of mesh points to use in the random model
268 //                      parameter selection process.  The upper limit is
269 //                      set by parameter nmax_integ in the accompanying
270 //                      include file 'Parameter_values.inc' which is presently
271 //                      set at 100, but the true upper limit for n_integ_pts
272 //                      is nmax_integ - 1 = 99. 
273 //   (18) n_scan_pts  - Number of mesh points to use to scan the (pt,y)
274 //                      dependence of the model distributions looking for
275 //                      the maximum value.  The 2-D grid has
276 //                      n_scan_pts * n_scan_pts points; no limit to size of
277 //                      n_scan_pts.
278 //   (19) irand       - Starting random number seed.
279 //
280 //**************************************************************************
281 //   FOR MODEL_TYPE = 1,2,3 or 4:
282 //   Input the following 7 lines for each particle type; repeat these
283 //   set of lines n_pid_type times:
284 //
285 //        (a) gpid - Geant Particle ID code number
286 //        (b) mult_mean,mult_variance_control - Mean multiplicity and
287 //                                              variance control where:
288 //            mult_variance_control = 0 for no variance in multiplicity 
289 //            mult_variance_control = 1 to allow Poisson distribution for
290 //                                      particle multiplicities for all events.
291 //            Note that a hard limit exists for the maximum possible
292 //            multiplicity for a given particle type per event.  This is
293 //            determined by parameter factorial_max in accompanying include
294 //            file 'common_facfac.inc' and is presently set at 10000.
295 //        (c) Temp_mean, Temp_stdev - Temperature parameter mean (in GeV)
296 //            and standard deviation (Gaussian distribution assumed).
297 //        (d) sigma_mean, sigma_stdev - Rapidity distribution width (sigma)
298 //            parameter mean and standard deviation (Gaussian distribution
299 //            assumed).
300 //        (e) expvel_mean, expvel_stdev - S. Pratt expansion velocity
301 //            (in units of c) mean and standard deviation (Gaussian 
302 //            distribution assumed).
303 //        (f) Vn_mean(k);  k=1,4  - Anisotropic flow parameters, mean values
304 //                                  for Fourier component n=1.
305 //        (g) Vn_stdev(k); k=1,4  - Anisotropic flow parameters, std.dev.
306 //                                  values for Fourier component n=1.
307 //
308 //            Repeat the last two lines of input for remaining Fourier
309 //            components n=2,3...6.  Include all 6 sets of parameters
310 //            even if these are not used by the model for Vn(pt,y) (set
311 //            unused parameter means and std.dev. to 0.0).  List 4 values
312 //            on every line, even though for n=even the 4th quantity is
313 //            not used.
314 //
315 //**************************************************************************
316 //   FOR MODEL_TYPE = 5 input the following set of lines for each particle
317 //                      type; repeat these n_pid_type times.
318 //
319 //        (a) gpid - Geant Particle ID code number
320 //        (b) mult_mean,mult_variance_control - Mean multiplicity and
321 //                                              variance control where:
322 //            mult_variance_control = 0 for no variance in multiplicity
323 //            mult_variance_control = 1 to allow Poisson distribution for
324 //                                      particle multiplicities for all events.
325 //        (c) pt_start, eta_start - minimum starting values for pt, eta 
326 //                                  input for the bin-by-bin distributions.
327 //        (d) n_pt_bins, n_eta_bins - # input pt and eta bins.
328 //        (e) delta_pt, pt_bin - pt bin size and function value, repeat for
329 //                               each pt bin.
330 //        (f) delta_eta, eta_bin - eta bin size and function value, repeat
331 //                                 for each eta bin.
332 //        (g) Vn_mean(k);  k=1,4  - Anisotropic flow parameters, mean values
333 //                                  for Fourier component n=1.
334 //        (h) Vn_stdev(k); k=1,4  - Anisotropic flow parameters, std.dev.
335 //                                  values for Fourier component n=1.
336 //
337 //            Repeat the last two lines of input for remaining Fourier
338 //            components n=2,3...6.  Include all 6 sets of parameters
339 //            even if these are not used by the model for Vn(pt,y) (set
340 //            unused parameter means and std.dev. to 0.0).  List 4 values
341 //            on every line, even though for n=even the 4th quantity is
342 //            not used.
343 //
344 //        NOTE: The pt, eta ranges must fully include the requested ranges
345 //              in input #4 and 5 above; else the code execution will stop.
346 //
347 //        Also, variable bin sizes are permitted for the input distributions.
348 //
349 //        Also, this input distribution is used for all events in the run;
350 //        no fluctuations in this "parent" distribution are allowed from 
351 //        event-to-event.
352 //
353 //**************************************************************************
354 //   FOR MODEL_TYPE = 6 input the following set of lines for each particle
355 //                      type; repeat these n_pid_type times.
356 //
357 //        (a) gpid - Geant Particle ID code number
358 //        (b) mult_mean,mult_variance_control - Mean multiplicity and
359 //                                              variance control where:
360 //            mult_variance_control = 0 for no variance in multiplicity
361 //            mult_variance_control = 1 to allow Poisson distribution for
362 //                                      particle multiplicities for all events.
363 //        (c) pt_start, eta_start - minimum starting values for pt, eta
364 //                                  input for the bin-by-bin distributions.
365 //        (d) n_pt_bins, n_eta_bins - # input pt and eta bins.
366 //        (e) delta_pt - pt bin size, repeat for each pt bin. 
367 //        (f) delta_eta - eta bin size, repeat for each eta bin.
368 //        (g) i,j,pt_eta_bin(i,j) - read pt (index = i) and eta (index = j)
369 //                                  bin numbers and bin value for full 2D space
370 //        (h) Vn_mean(k);  k=1,4  - Anisotropic flow parameters, mean values
371 //                                  for Fourier component n=1.
372 //        (i) Vn_stdev(k); k=1,4  - Anisotropic flow parameters, std.dev.
373 //                                  values for Fourier component n=1.
374 //
375 //            Repeat the last two lines of input for remaining Fourier
376 //            components n=2,3...6.  Include all 6 sets of parameters
377 //            even if these are not used by the model for Vn(pt,y) (set
378 //            unused parameter means and std.dev. to 0.0).  List 4 values
379 //            on every line, even though for n=even the 4th quantity is
380 //            not used.
381 //
382 //        NOTE: The pt, eta ranges must fully include the requested ranges
383 //              in input #4 and 5 above; else the code execution will stop.
384 //
385 //        Also, variable bin sizes are permitted for the input distributions.
386 //
387 //        Also, this input distribution is used for all events in the run;
388 //        no fluctuations in this "parent" distribution are allowed from
389 //        event-to-event.
390 //
391 ///////////////////////////////////////////////////////////////////////////////
392
393
394 #include <Riostream.h>
395 #include <TParticle.h>
396 #include <TClonesArray.h>
397
398 #include "TMevSim.h"
399 #define IN_TMEVSIM_CXX
400 #include "TMevSimPartTypeParams.h"
401 #include "TParticle.h"
402
403 #ifndef WIN32
404 # define multgen multgen_
405 # define type_of_call
406 #else
407 # define multgen MULTGEN
408 # define type_of_call _stdcall
409 #endif
410
411
412 ClassImp(TMevSim)
413
414
415 extern "C" void type_of_call multgen();
416
417 //______________________________________________________________________________
418 TMevSim::TMevSim(Int_t nEvents, Int_t modelType, Int_t reacPlaneCntrl,
419            Float_t psiRMean, Float_t psiRStDev, Float_t multFacMean, Float_t multFacStDev,
420            Float_t ptCutMin, Float_t ptCutMax, Float_t etaCutMin, Float_t etaCutMax, 
421            Float_t phiCutMin, Float_t phiCutMax, Int_t irand) : TGenerator("MevSim", "MevSim")
422 {
423 // TMevSim constructor: initializes all the event-wide variables of MevSim with
424 // user supplied values, or with the default ones (declared in the header file).
425 // It also allocates space for the array which will store parameters specific to
426 // each particle species.
427 // Caution: Setting nEvents > 1 will have no effect, since only the last generated 
428 // event will be stored in POUT COMMON, and therefore only one event can be 
429 // accessible at a time. 
430  
431    fNEvents = nEvents;
432    fModelType = modelType;
433    fReacPlaneCntrl = reacPlaneCntrl;
434    fPsiRMean = psiRMean;
435    fPsiRStDev = psiRStDev;
436    fMultFacMean = multFacMean;
437    fMultFacStDev = multFacStDev;
438    fPtCutMin = ptCutMin;
439    fPtCutMax = ptCutMax;
440    fEtaCutMin = etaCutMin;
441    fEtaCutMax = etaCutMax;
442    fPhiCutMin = phiCutMin;
443    fPhiCutMax = phiCutMax;
444    fNStDevMult = fNStDevTemp = fNStDevSigma = fNStDevExpVel = fNStdDevPSIr = fNStDevVn = fNStDevMultFac = 3.0;
445    fNIntegPts = 100;
446    fNScanPts = 100;
447    firand = irand;
448    fParticleTypeParameters = new TClonesArray("TMevSimPartTypeParams",10);
449    fNPDGCodes = 0;
450    DefineParticles();
451 }
452 //______________________________________________________________________________
453 TMevSim::~TMevSim()
454 {
455 // TMevSim destructor: destroys the object and all the particle information stored
456 // in the list.
457    
458   if (fParticleTypeParameters) {
459     fParticleTypeParameters->Clear();
460     delete fParticleTypeParameters;
461     fParticleTypeParameters = 0;
462   }
463 }
464 //______________________________________________________________________________
465 TMevSim::TMevSim(TMevSim& mevsim) : TGenerator(mevsim) {
466 // The copy constructor
467    
468    *this = mevsim;
469 }
470 //______________________________________________________________________________
471
472 TMevSim& TMevSim::operator=(const TMevSim& mevsim) {
473 // An assignment operator: initializes all the event-wide variables of MevSim with
474 // the ones from a copied object. It also copies the parameters specific to
475 // each particle species.
476    
477    fNEvents = mevsim.GetNEvents();
478    fModelType = mevsim.GetModelType();
479    fReacPlaneCntrl = mevsim.GetReacPlaneCntrl();
480    fPsiRMean = mevsim.GetPsiRMean();
481    fPsiRStDev = mevsim.GetPsiRStDev();
482    fMultFacMean = mevsim.GetMultFacMean();
483    fMultFacStDev = mevsim.GetMultFacStDev();
484    fPtCutMin = mevsim.GetPtCutMin();
485    fPtCutMax = mevsim.GetPtCutMax();
486    fEtaCutMin = mevsim.GetEtaCutMin();
487    fEtaCutMax = mevsim.GetEtaCutMax();
488    fPhiCutMin = mevsim.GetPhiCutMin();
489    fPhiCutMax = mevsim.GetPhiCutMax();
490    fNStDevMult = mevsim.GetNStDevMult();
491    fNStDevTemp = mevsim.GetNStDevTemp();
492    fNStDevSigma =GetNStDevSigma();
493    fNStDevExpVel = mevsim.GetNStDevExpVel();
494    fNStdDevPSIr = mevsim.GetNStDevPSIr(); 
495    fNStDevVn =  mevsim.GetNStDevVn();
496    fNStDevMultFac =  mevsim.GetNStDevMultFac();
497    fNIntegPts = mevsim.GetNintegPts();
498    fNScanPts = mevsim.GetNScanPts();
499    firand = mevsim.firand;
500    fParticleTypeParameters = new TClonesArray("TMevSimPartTypeParams",mevsim.GetNPidTypes());
501    for (int i=0; i< mevsim.GetNPidTypes(); i++) 
502      {  
503         TMevSimPartTypeParams *temp = 0;
504         mevsim.GetPartTypeParamsByIndex(i,temp);
505         fParticleTypeParameters->AddLast(temp);
506      }
507    DefineParticles();
508    return (*this);
509 }
510 //______________________________________________________________________________
511 void        TMevSim::Initialize() {
512 // TMevSim initialization: creates an input file for the FORTRAN
513 // program MevSim. Converts all the event-wide information and particle
514 // specific information to the format readable by MevSim and writes it
515 // to disk in current directory.
516 // Caution: At least one TMevSimPartTypeParams object must be created and 
517 // added to the collection before event generation can start.
518    
519   TMevSimPartTypeParams * params = 0;
520   
521
522    ofstream *file = new ofstream("mult_gen.in",ios::out | ios::trunc);
523    // Write out the parameters to the pramameter file
524    *file << "  " << fNEvents << "  ! Number of Events \n";
525    *file << "  " << GetNPidTypes() << " \n";
526    *file << "  " << fModelType << " \n";
527    *file << "  " << fReacPlaneCntrl << " \n";
528    file->setf(ios::showpoint);
529    *file << "  " << fPsiRMean << "   " << fPsiRStDev << " \n";
530    *file << "  " << fMultFacMean << "   " << fMultFacStDev << " \n";
531    *file << "  " << fPtCutMin << "   " << fPtCutMax << " \n";
532    *file << "  " << fEtaCutMin << "   " << fEtaCutMax << " \n";
533    *file << "  " << fPhiCutMin << "   " << fPhiCutMax << " \n";
534    *file << "  " << fNStDevMult << " \n";
535    *file << "  " << fNStDevTemp << " \n";
536    *file << "  " << fNStDevSigma << " \n";
537    *file << "  " << fNStDevExpVel << " \n";
538    *file << "  " << fNStdDevPSIr << " \n";
539    *file << "  " << fNStDevVn << " \n";
540    *file << "  " << fNStDevMultFac << " \n";
541    *file << "  " << fNIntegPts << " \n";
542    *file << "  " << fNScanPts << " \n";
543    *file << "  " << firand << " \n";
544    // Write out particle specific information
545    for (Int_t i=0; i< (fParticleTypeParameters->GetLast() + 1); i++) {
546
547      params = (TMevSimPartTypeParams *) ((*fParticleTypeParameters)[i]);
548      
549      *file << "  " << params->GetGPid() << "       ! Particle GEANT Pid \n";
550      *file << "  " << params->GetMultMean() << "  " << params->GetMultVarianceControl() << "  \n";
551      *file << "  " << params->GetTempMean() << "  " << params->GetTempStDev() << "  \n";
552      *file << "  " << params->GetSigmaMean() << "  " << params->GetSigmaStDev() << "  \n";
553      *file << "  " << params->GetExpVelMean() << "  " << params->GetExpVelStDev() << "  \n";
554
555      for (Int_t cnt1 = 0; cnt1 < NFLOWTERMS; cnt1++) {
556        *file << "  ";
557        Int_t cnt2;
558        for (cnt2 = 0; cnt2 < 4; cnt2++) *file << params->GetVnMeanComponent(cnt1, cnt2) << "  ";
559        *file << "  \n  ";
560        for (cnt2 = 0; cnt2 < 4; cnt2++) *file << params->GetVnStDevComponent(cnt1, cnt2) << "  ";
561        *file << "  \n";
562      }
563    }
564    file->close();
565
566 }
567 //______________________________________________________________________________
568 void TMevSim::GenerateEvent() {
569 // Generates one MevSim event. TMevSim::Initialize() must be called prior
570 // to calling this function.
571    
572    Info("GenerateEvent","Calling FORTRAN multgen()");
573    multgen();
574 }
575
576 //______________________________________________________________________________
577 Int_t TMevSim::ImportParticles(TClonesArray *particles, Option_t */*option*/)
578 {
579 // Read in particles created by MevSim into the TClonesArray(). The Initialize()
580 // and GenrateEvent() functions must be called prior to calling this funtion.
581 // The particles are read from the COMMON POUT. Right now the only provided 
582 // information is Geant PID, 3 momentum components and the energy of the particle.
583    
584    if (particles == 0) return 0;
585    TClonesArray &aParticles = *particles;
586    aParticles.Clear();
587    
588    Int_t totpart = 0;
589    for (Int_t nrpidtype=0; nrpidtype < (fParticleTypeParameters->GetLast() + 1); nrpidtype++) {
590       Int_t nrpart = 0;
591       Int_t pidcode = ((TMevSimPartTypeParams *) (*fParticleTypeParameters)[nrpidtype])->GetGPid();
592       while ((TRACK.pout[(4*nrpart+3)*NPID+nrpidtype] > 0.0) || (TRACK.pout[(4*nrpart)*NPID+nrpidtype] != 0.0)) {
593          int poffset = 4*nrpart*NPID+nrpidtype;
594          Float_t px = TRACK.pout[poffset];
595          poffset += NPID;
596          Float_t py = TRACK.pout[poffset];
597          poffset += NPID;
598          Float_t pz = TRACK.pout[poffset];
599          poffset += NPID;
600          Float_t mass = TRACK.pout[poffset];
601          new(aParticles[totpart+nrpart]) TParticle(
602                                           PDGFromId(pidcode),  // Get the PDG ID from GEANT ID
603                                           0,
604                                           0,
605                                           0,
606                                           0,
607                                           0,
608                                           px,
609                                           py,
610                                           pz,
611                                           sqrt(mass*mass+px*px+py*py+pz*pz),                               
612                                           0,
613                                           0,
614                                           0,
615                                           0);
616          nrpart++;
617       }
618       totpart += nrpart;
619    }
620    return totpart;
621 }
622 //______________________________________________________________________________
623 TObjArray * TMevSim::ImportParticles(Option_t */*option*/)
624 {
625 // Read in particles created by MevSim into the TClonesArray(). The Initialize()
626 // and GenrateEvent() functions must be called prior to calling this funtion.
627 // The particles are read from the COMMON POUT. Right now the only provided 
628 // information is Geant PID, 3 momentum components and the energy of the particle.
629    
630    fParticles->Clear();
631    
632    for (Int_t nrpidtype=0; nrpidtype < (fParticleTypeParameters->GetLast() + 1); nrpidtype++) {
633       Int_t nrpart = 0;
634       Int_t pidcode = ((TMevSimPartTypeParams *) (*fParticleTypeParameters)[nrpidtype])->GetGPid();
635       while ((TRACK.pout[(4*nrpart+3)*NPID+nrpidtype] > 0.0) || (TRACK.pout[(4*nrpart)*NPID+nrpidtype] != 0.0)) {
636          int poffset = 4*nrpart*NPID+nrpidtype;
637          Float_t px = TRACK.pout[poffset];
638          poffset += NPID;
639          Float_t py = TRACK.pout[poffset];
640          poffset += NPID;
641          Float_t pz = TRACK.pout[poffset];
642          poffset += NPID;
643          Float_t mass = TRACK.pout[poffset];
644          TParticle * p = new TParticle(
645                                           PDGFromId(pidcode),  // Get the PDG ID from GEANT ID
646                                           0,
647                                           0,
648                                           0,
649                                           0,
650                                           0,
651                                           px,
652                                           py,
653                                           pz,
654                                           sqrt(mass*mass+px*px+py*py+pz*pz),                               
655                                           0,
656                                           0,
657                                           0,
658                                           0);
659          fParticles->Add(p);
660          nrpart++;
661       }
662    }
663    return fParticles;
664 }
665 //______________________________________________________________________________
666 void        TMevSim::SetNEvents(Int_t nEvents ) { 
667 // Sets the number of events to be generated by MevSim.
668 // Caution: Setting nEvents > 1 will have no effect, since only the last generated 
669 // event will be stored in POUT COMMON, and therefore only one event can be 
670 // accessible at a time. 
671
672    fNEvents = nEvents;
673 }
674 //______________________________________________________________________________
675 Int_t       TMevSim::GetNEvents() const {
676    return fNEvents;
677 }
678 //______________________________________________________________________________
679 Int_t       TMevSim::GetNPidTypes() const {
680    return fParticleTypeParameters->GetLast()+1;
681 }
682 //______________________________________________________________________________
683 void        TMevSim::SetModelType(Int_t modelType) {
684    fModelType  = modelType;
685 }
686 //______________________________________________________________________________
687 Int_t       TMevSim::GetModelType() const {
688    return fModelType;
689 }
690 //______________________________________________________________________________
691 void        TMevSim::SetReacPlaneCntrl(Int_t reacPlaneCntrl) {
692    fReacPlaneCntrl = reacPlaneCntrl; 
693 }
694 //______________________________________________________________________________
695 Int_t       TMevSim::GetReacPlaneCntrl() const {
696    return fReacPlaneCntrl;
697 }
698 //______________________________________________________________________________
699 void        TMevSim::SetPsiRParams(Float_t psiRMean,  Float_t psiRStDev) {
700    fPsiRMean = psiRMean;
701    fPsiRStDev = psiRStDev;
702 }
703 //______________________________________________________________________________
704 Float_t     TMevSim::GetPsiRMean() const { 
705    return fPsiRMean;
706 }
707 //______________________________________________________________________________
708 Float_t     TMevSim::GetPsiRStDev() const { 
709    return fPsiRStDev;
710 }
711 //______________________________________________________________________________
712 void        TMevSim::SetMultFacParams(Float_t multFacMean,  Float_t multFacStDev) {
713    fMultFacMean = multFacMean;
714    fMultFacStDev = multFacStDev;
715 }
716 //______________________________________________________________________________
717 Float_t     TMevSim::GetMultFacMean() const { 
718    return fMultFacMean;     
719 }
720 //______________________________________________________________________________
721 Float_t     TMevSim::GetMultFacStDev() const { 
722    return fMultFacStDev;
723 }
724 //______________________________________________________________________________
725 void        TMevSim::SetPtCutRange(Float_t ptCutMin,  Float_t ptCutMax) { 
726    fPtCutMin = ptCutMin;
727    fPtCutMax = ptCutMax;
728 }
729 //______________________________________________________________________________
730 Float_t     TMevSim::GetPtCutMin() const { 
731    return fPtCutMin;
732 }
733 //______________________________________________________________________________
734 Float_t     TMevSim::GetPtCutMax() const { 
735    return fPtCutMax;
736 }
737 //______________________________________________________________________________
738 void        TMevSim::SetEtaCutRange(Float_t etaCutMin,  Float_t etaCutMax) {     fEtaCutMin = etaCutMin;
739    fEtaCutMax = etaCutMax;
740 }
741
742 //______________________________________________________________________________
743    Float_t     TMevSim::GetEtaCutMin() const { 
744      return fEtaCutMin;
745 }
746 //______________________________________________________________________________
747    Float_t     TMevSim::GetEtaCutMax() const {
748      return fEtaCutMax;
749 }
750 //______________________________________________________________________________
751 void        TMevSim::SetPhiCutRange(Float_t phiCutMin,  Float_t phiCutMax) {
752    fPhiCutMin = phiCutMin;
753    fPhiCutMax = phiCutMax;
754 }
755 //______________________________________________________________________________
756 Float_t     TMevSim::GetPhiCutMin() const { 
757    return fPhiCutMin;
758 }
759 //______________________________________________________________________________
760 Float_t     TMevSim::GetPhiCutMax() const { 
761    return fPhiCutMax;
762 }
763 //______________________________________________________________________________
764 void        TMevSim::SetNStDevMult(Float_t nStDevMult) { 
765    fNStDevMult = nStDevMult;
766 }
767 //______________________________________________________________________________
768 Float_t       TMevSim::GetNStDevMult() const { 
769    return  fNStDevMult;
770 }
771 //______________________________________________________________________________
772 void        TMevSim::SetNStDevTemp(Float_t nStDevTemp) { 
773    fNStDevTemp = nStDevTemp;
774 }
775 //______________________________________________________________________________
776 Float_t       TMevSim::GetNStDevTemp() const { 
777    return fNStDevTemp;
778 }
779 //______________________________________________________________________________
780 void        TMevSim::SetNStDevSigma(Float_t nStDevSigma) { 
781    fNStDevSigma = nStDevSigma;
782 }
783 //______________________________________________________________________________
784 Float_t       TMevSim::GetNStDevSigma() const { 
785    return fNStDevSigma;  
786 }
787 //______________________________________________________________________________
788 void        TMevSim::SetNStDevExpVel(Float_t nStDevExpVel) { 
789    fNStDevExpVel = nStDevExpVel;
790 }
791 //______________________________________________________________________________
792 Float_t       TMevSim::GetNStDevExpVel() const { 
793    return fNStDevExpVel;
794 }
795 //______________________________________________________________________________
796 void        TMevSim::SetNStDevPSIr(Float_t nStDevPSIr) { 
797    fNStdDevPSIr = nStDevPSIr;
798 }
799 //______________________________________________________________________________
800 Float_t       TMevSim::GetNStDevPSIr() const { 
801    return fNStdDevPSIr;
802 }
803 //______________________________________________________________________________
804 void        TMevSim::SetNStDevVn(Float_t nStDevVn) { 
805    fNStDevVn = nStDevVn;
806 }
807 //______________________________________________________________________________
808 Float_t       TMevSim::GetNStDevVn() const { 
809    return fNStDevVn; 
810 }
811 //______________________________________________________________________________
812 void        TMevSim::SetNStDevMultFac(Float_t nStDevMultFac) { 
813    fNStDevMultFac = nStDevMultFac;
814 }
815 //______________________________________________________________________________
816 Float_t       TMevSim::GetNStDevMultFac() const { 
817    return fNStDevMultFac;  
818 }
819 //______________________________________________________________________________
820 void        TMevSim::SetNIntegPts(Int_t nIntegPts) { 
821    fNIntegPts = nIntegPts;
822 }
823 //______________________________________________________________________________
824 Int_t       TMevSim::GetNintegPts() const { 
825    return fNIntegPts;
826 }
827 //______________________________________________________________________________
828 void        TMevSim::SetNScanPts(Int_t nScanPts) { 
829    fNScanPts = nScanPts;
830 }
831 //______________________________________________________________________________
832 Int_t       TMevSim::GetNScanPts() const { 
833    return fNScanPts;  
834 }
835 //______________________________________________________________________________
836 void        TMevSim::AddPartTypeParams(TMevSimPartTypeParams *params) {
837 // Add the particle specied parameters and the end of the list.
838    
839   //cout << params << " " <<  fParticleTypeParameters <<  endl;
840
841   //fParticleTypeParameters->Dump(); 
842   params->Dump();
843   
844   Int_t last = fParticleTypeParameters->GetLast();
845   new ((*fParticleTypeParameters)[last+1]) TMevSimPartTypeParams(*params);
846 }
847 //______________________________________________________________________________
848 void        TMevSim::SetPartTypeParams(Int_t index, TMevSimPartTypeParams *params) 
849 {
850 // Create the new copy particle species parameters provided by params, and store
851 // them in the position designated by index. 
852    
853    *((TMevSimPartTypeParams *) ((*fParticleTypeParameters)[index])) = *params;
854 }
855 //______________________________________________________________________________
856 void TMevSim::GetPartTypeParamsByIndex(Int_t index, TMevSimPartTypeParams *params) const
857 {
858 // Return the particle parameters stored in the list at the postion index.   
859 // Returns NULL if index is out of bounds.
860    
861    if ((index < fParticleTypeParameters->GetLast()) && (index >= 0))
862         params = (TMevSimPartTypeParams *) (*fParticleTypeParameters)[index];
863    else
864      params = NULL;
865 }
866 //______________________________________________________________________________
867 void TMevSim::GetPartTypeParamsByGPid(Int_t gpid, TMevSimPartTypeParams *params) const
868 {
869 // Return the particle parameters for the particle with Geant PID gpid.
870 // Returns NULL if the parameters for such particle do not exist in the list.   
871    
872    Int_t i = -1;
873    
874    while (++i <= fParticleTypeParameters->GetLast())
875      {
876         if (((TMevSimPartTypeParams *) (*fParticleTypeParameters)[i])->GetGPid() == gpid)
877           {
878              params = (TMevSimPartTypeParams *) (*fParticleTypeParameters)[i];
879              return;
880           }
881      }
882    params = NULL;
883    return;
884 }
885 //_____________________________________________________________________________
886 Int_t TMevSim::PDGFromId(Int_t id) const 
887 {
888   //
889   // Return PDG code and pseudo ENDF code from Geant3 code
890   //
891   if(id>0 && id<fNPDGCodes) return fPDGCode[id];
892   else return -1;
893 }
894 //_____________________________________________________________________________
895 void TMevSim::DefineParticles() 
896 {
897   //
898   // Load standard numbers for GEANT particles and PDG conversion
899   fPDGCode[fNPDGCodes++]=-99;   //  0 = unused location
900   fPDGCode[fNPDGCodes++]=22;    //  1 = photon
901   fPDGCode[fNPDGCodes++]=-11;   //  2 = positron
902   fPDGCode[fNPDGCodes++]=11;    //  3 = electron
903   fPDGCode[fNPDGCodes++]=12;    //  4 = neutrino e
904   fPDGCode[fNPDGCodes++]=-13;   //  5 = muon +
905   fPDGCode[fNPDGCodes++]=13;    //  6 = muon -
906   fPDGCode[fNPDGCodes++]=111;   //  7 = pi0
907   fPDGCode[fNPDGCodes++]=211;   //  8 = pi+
908   fPDGCode[fNPDGCodes++]=-211;  //  9 = pi-
909   fPDGCode[fNPDGCodes++]=130;   // 10 = Kaon Long
910   fPDGCode[fNPDGCodes++]=321;   // 11 = Kaon +
911   fPDGCode[fNPDGCodes++]=-321;  // 12 = Kaon -
912   fPDGCode[fNPDGCodes++]=2112;  // 13 = Neutron
913   fPDGCode[fNPDGCodes++]=2212;  // 14 = Proton
914   fPDGCode[fNPDGCodes++]=-2212; // 15 = Anti Proton
915   fPDGCode[fNPDGCodes++]=310;   // 16 = Kaon Short
916   fPDGCode[fNPDGCodes++]=221;   // 17 = Eta
917   fPDGCode[fNPDGCodes++]=3122;  // 18 = Lambda
918   fPDGCode[fNPDGCodes++]=3222;  // 19 = Sigma +
919   fPDGCode[fNPDGCodes++]=3212;  // 20 = Sigma 0
920   fPDGCode[fNPDGCodes++]=3112;  // 21 = Sigma -
921   fPDGCode[fNPDGCodes++]=3322;  // 22 = Xi0
922   fPDGCode[fNPDGCodes++]=3312;  // 23 = Xi-
923   fPDGCode[fNPDGCodes++]=3334;  // 24 = Omega-
924   fPDGCode[fNPDGCodes++]=-2112; // 25 = Anti Proton
925   fPDGCode[fNPDGCodes++]=-3122; // 26 = Anti Proton
926   fPDGCode[fNPDGCodes++]=-3222; // 27 = Anti Sigma -
927   fPDGCode[fNPDGCodes++]=-3212; // 28 = Anti Sigma 0
928   fPDGCode[fNPDGCodes++]=-3112; // 29 = Anti Sigma 0
929   fPDGCode[fNPDGCodes++]=-3322; // 30 = Anti Xi 0
930   fPDGCode[fNPDGCodes++]=-3312; // 31 = Anti Xi +
931   fPDGCode[fNPDGCodes++]=-3334; // 32 = Anti Omega +
932 }
933