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