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