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