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