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