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