1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 Revision 1.2 2001/03/28 07:32:51 hristov
19 Loop variables declared only once, old style include (HP,Sun)
21 Revision 1.1 2001/03/25 10:15:23 morsch
22 Root interface to MevSim code as TGenerator realisation (Sylwester Radomski et al.)
26 ////////////////////////////////////////////////////////////////////////////
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.
35 // For The STAR Collaboration
39 // The University of Texas at Austin
40 // Austin, Texas 78712
42 // ray@physics.utexas.edu
44 // Ron Longacre email:
47 ////////////////////////////////////////////////////////////////////////////
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.
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,
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)]}]
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.
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.
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
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
103 // (4) assume uniformly distributed, random values for the reaction
104 // plane angles from 0 to 360 deg., for each event in the run.
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.
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.
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
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.
157 // An interpolation subroutine from Rubin Landau, Oregon State Univ.,
158 // is used to do this interpolation; it involves uneven mesh point
161 // The method for generating the particle momenta uses the
162 // standard random elimination method and involves the following
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
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
186 // For model_type = 5,6 (see following) which are input bin-by-bin
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
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.
214 // III. DESCRIPTION OF THE INPUT:
217 // The input is described below in the 'read' statements and also in
218 // the sample input file. Some additional comments are as follows:
220 // (1) n_events - Selected number of events in run. Can be anything
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
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
286 // (19) irand - Starting random number seed.
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:
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
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.
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
323 //**************************************************************************
324 // FOR MODEL_TYPE = 5 input the following set of lines for each particle
325 // type; repeat these n_pid_type times.
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
338 // (f) delta_eta, eta_bin - eta bin size and function value, repeat
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.
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
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.
355 // Also, variable bin sizes are permitted for the input distributions.
357 // Also, this input distribution is used for all events in the run;
358 // no fluctuations in this "parent" distribution are allowed from
361 //**************************************************************************
362 // FOR MODEL_TYPE = 6 input the following set of lines for each particle
363 // type; repeat these n_pid_type times.
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.
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
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.
393 // Also, variable bin sizes are permitted for the input distributions.
395 // Also, this input distribution is used for all events in the run;
396 // no fluctuations in this "parent" distribution are allowed from
399 ///////////////////////////////////////////////////////////////////////////////
408 #include "MevSimCommon.h"
409 #include "TParticle.h"
413 # define multgen multgen_
414 # define type_of_call
416 # define multgen MULTGEN
417 # define type_of_call _stdcall
424 extern "C" void type_of_call multgen();
426 //______________________________________________________________________________
427 TMevSim::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")
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.
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;
457 fParticleTypeParameters = new TClonesArray("TMevSimPartTypeParams",10);
461 //______________________________________________________________________________
464 // TMevSim destructor: destroys the object and all the particle information stored
467 if (fParticleTypeParameters) {
468 fParticleTypeParameters->Clear();
469 delete fParticleTypeParameters;
470 fParticleTypeParameters = 0;
473 //______________________________________________________________________________
474 TMevSim::TMevSim(TMevSim& mevsim) {
475 // The copy constructor
479 //______________________________________________________________________________
480 TMevSim& 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.
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++)
511 TMevSimPartTypeParams *temp = 0;
512 mevsim.GetPartTypeParamsByIndex(i,temp);
513 fParticleTypeParameters->AddLast(temp);
518 //______________________________________________________________________________
519 void 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.
527 TMevSimPartTypeParams * params = 0;
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++) {
555 params = (TMevSimPartTypeParams *) ((*fParticleTypeParameters)[i]);
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";
563 for (Int_t cnt1 = 0; cnt1 < NFLOWTERMS; cnt1++) {
566 for (cnt2 = 0; cnt2 < 4; cnt2++) *file << params->GetVnMeanComponent(cnt1, cnt2) << " ";
568 for (cnt2 = 0; cnt2 < 4; cnt2++) *file << params->GetVnStDevComponent(cnt1, cnt2) << " ";
575 //______________________________________________________________________________
576 void TMevSim::GenerateEvent() {
577 // Generates one MevSim event. TMevSim::Initialize() must be called prior
578 // to calling this function.
580 cout << "Calling FORTRAN multgen()" << endl;
584 //______________________________________________________________________________
585 Int_t TMevSim::ImportParticles(TClonesArray *particles, Option_t *option)
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.
592 if (particles == 0) return 0;
593 TClonesArray &Particles = *particles;
597 for (Int_t nrpidtype=0; nrpidtype < (fParticleTypeParameters->GetLast() + 1); nrpidtype++) {
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];
604 Float_t py = TRACK.pout[poffset];
606 Float_t pz = TRACK.pout[poffset];
608 Float_t mass = TRACK.pout[poffset];
609 new(Particles[totpart+nrpart]) TParticle(
610 PDGFromId(pidcode), // Get the PDG ID from GEANT ID
619 sqrt(mass*mass+px*px+py*py+pz*pz),
630 //______________________________________________________________________________
631 void 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.
639 //______________________________________________________________________________
640 Int_t TMevSim::GetNEvents() const {
643 //______________________________________________________________________________
644 Int_t TMevSim::GetNPidTypes() const {
645 return fParticleTypeParameters->GetLast()+1;
647 //______________________________________________________________________________
648 void TMevSim::SetModelType(Int_t modelType) {
649 fModelType = modelType;
651 //______________________________________________________________________________
652 Int_t TMevSim::GetModelType() const {
655 //______________________________________________________________________________
656 void TMevSim::SetReacPlaneCntrl(Int_t reacPlaneCntrl) {
657 fReacPlaneCntrl = reacPlaneCntrl;
659 //______________________________________________________________________________
660 Int_t TMevSim::GetReacPlaneCntrl() const {
661 return fReacPlaneCntrl;
663 //______________________________________________________________________________
664 void TMevSim::SetPsiRParams(Float_t psiRMean, Float_t psiRStDev) {
665 fPsiRMean = psiRMean;
666 fPsiRStDev = psiRStDev;
668 //______________________________________________________________________________
669 Float_t TMevSim::GetPsiRMean() const {
672 //______________________________________________________________________________
673 Float_t TMevSim::GetPsiRStDev() const {
676 //______________________________________________________________________________
677 void TMevSim::SetMultFacParams(Float_t multFacMean, Float_t multFacStDev) {
678 fMultFacMean = multFacMean;
679 fMultFacStDev = multFacStDev;
681 //______________________________________________________________________________
682 Float_t TMevSim::GetMultFacMean() const {
685 //______________________________________________________________________________
686 Float_t TMevSim::GetMultFacStDev() const {
687 return fMultFacStDev;
689 //______________________________________________________________________________
690 void TMevSim::SetPtCutRange(Float_t ptCutMin, Float_t ptCutMax) {
691 fPtCutMin = ptCutMin;
692 fPtCutMax = ptCutMax;
694 //______________________________________________________________________________
695 Float_t TMevSim::GetPtCutMin() const {
698 //______________________________________________________________________________
699 Float_t TMevSim::GetPtCutMax() const {
702 //______________________________________________________________________________
703 void TMevSim::SetEtaCutRange(Float_t etaCutMin, Float_t etaCutMax) { fEtaCutMin = etaCutMin;
704 fEtaCutMax = etaCutMax;
707 //______________________________________________________________________________
708 Float_t TMevSim::GetEtaCutMin() const {
711 //______________________________________________________________________________
712 Float_t TMevSim::GetEtaCutMax() const {
715 //______________________________________________________________________________
716 void TMevSim::SetPhiCutRange(Float_t phiCutMin, Float_t phiCutMax) {
717 fPhiCutMin = phiCutMin;
718 fPhiCutMax = phiCutMax;
720 //______________________________________________________________________________
721 Float_t TMevSim::GetPhiCutMin() const {
724 //______________________________________________________________________________
725 Float_t TMevSim::GetPhiCutMax() const {
728 //______________________________________________________________________________
729 void TMevSim::SetNStDevMult(Float_t nStDevMult) {
730 fNStDevMult = nStDevMult;
732 //______________________________________________________________________________
733 Float_t TMevSim::GetNStDevMult() const {
736 //______________________________________________________________________________
737 void TMevSim::SetNStDevTemp(Float_t nStDevTemp) {
738 fNStDevTemp = nStDevTemp;
740 //______________________________________________________________________________
741 Float_t TMevSim::GetNStDevTemp() const {
744 //______________________________________________________________________________
745 void TMevSim::SetNStDevSigma(Float_t nStDevSigma) {
746 fNStDevSigma = nStDevSigma;
748 //______________________________________________________________________________
749 Float_t TMevSim::GetNStDevSigma() const {
752 //______________________________________________________________________________
753 void TMevSim::SetNStDevExpVel(Float_t nStDevExpVel) {
754 fNStDevExpVel = nStDevExpVel;
756 //______________________________________________________________________________
757 Float_t TMevSim::GetNStDevExpVel() const {
758 return fNStDevExpVel;
760 //______________________________________________________________________________
761 void TMevSim::SetNStDevPSIr(Float_t nStDevPSIr) {
762 fNStdDevPSIr = nStDevPSIr;
764 //______________________________________________________________________________
765 Float_t TMevSim::GetNStDevPSIr() const {
768 //______________________________________________________________________________
769 void TMevSim::SetNStDevVn(Float_t nStDevVn) {
770 fNStDevVn = nStDevVn;
772 //______________________________________________________________________________
773 Float_t TMevSim::GetNStDevVn() const {
776 //______________________________________________________________________________
777 void TMevSim::SetNStDevMultFac(Float_t nStDevMultFac) {
778 fNStDevMultFac = nStDevMultFac;
780 //______________________________________________________________________________
781 Float_t TMevSim::GetNStDevMultFac() const {
782 return fNStDevMultFac;
784 //______________________________________________________________________________
785 void TMevSim::SetNIntegPts(Int_t nIntegPts) {
786 fNIntegPts = nIntegPts;
788 //______________________________________________________________________________
789 Int_t TMevSim::GetNintegPts() const {
792 //______________________________________________________________________________
793 void TMevSim::SetNScanPts(Int_t nScanPts) {
794 fNScanPts = nScanPts;
796 //______________________________________________________________________________
797 Int_t TMevSim::GetNScanPts() const {
800 //______________________________________________________________________________
801 void TMevSim::AddPartTypeParams(TMevSimPartTypeParams *params) {
802 // Add the particle specied parameters and the end of the list.
804 //cout << params << " " << fParticleTypeParameters << endl;
806 //fParticleTypeParameters->Dump();
809 Int_t last = fParticleTypeParameters->GetLast();
810 new ((*fParticleTypeParameters)[last+1]) TMevSimPartTypeParams(*params);
812 //______________________________________________________________________________
813 void TMevSim::SetPartTypeParams(Int_t index, TMevSimPartTypeParams *params)
815 // Create the new copy particle species parameters provided by params, and store
816 // them in the position designated by index.
818 *((TMevSimPartTypeParams *) ((*fParticleTypeParameters)[index])) = *params;
820 //______________________________________________________________________________
821 void TMevSim::GetPartTypeParamsByIndex(Int_t index, TMevSimPartTypeParams *params)
823 // Return the particle parameters stored in the list at the postion index.
824 // Returns NULL if index is out of bounds.
826 if ((index < fParticleTypeParameters->GetLast()) && (index >= 0))
827 params = (TMevSimPartTypeParams *) (*fParticleTypeParameters)[index];
831 //______________________________________________________________________________
832 void TMevSim::GetPartTypeParamsByGPid(Int_t gpid, TMevSimPartTypeParams *params)
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.
839 while (++i <= fParticleTypeParameters->GetLast())
841 if (((TMevSimPartTypeParams *) (*fParticleTypeParameters)[i])->GetGPid() == gpid)
843 params = (TMevSimPartTypeParams *) (*fParticleTypeParameters)[i];
850 //_____________________________________________________________________________
851 Int_t TMevSim::PDGFromId(Int_t id) const
854 // Return PDG code and pseudo ENDF code from Geant3 code
856 if(id>0 && id<fNPDGCodes) return fPDGCode[id];
859 //_____________________________________________________________________________
860 void TMevSim::DefineParticles()
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 +