Pluto v6.3.0
Loading...
Searching...
No Matches
PHUrReader.cc
Go to the documentation of this file.
1
2// Read in urqmd f15 file (collision history) and put the particles
3// on the PLUTO bulk stack. Either dileptons (D0,D+,w,rho0,phi,eta,eta'
4// and pi0) can be decayed via shining (translated for original UrQMD
5// fortran routines, authors: Katharina Schmidt, Sascha Vogel,
6// Christian Sturm 2007-2009) into e+/e- or by PLUTO bulk decay.
7// Optional stable non dileptons can be written out. PHUrReader adds
8// a TClonesArray "Addon" to the PLUTO Tree that contains original
9// UrQMD infos about creation and absorbtion times, densities etc.
10//
11// USAGE: see urqmd_f15_input.C
13
14#include "PHUrReader.h"
15#include "TSystem.h"
16#include "PStaticData.h"
17#include "PDynamicData.h"
18
19#include "PDistributionManager.h"
20
21ClassImp(PHUrReader)
22
23
24TClonesArray *PHUrReader::fAdd = 0;
25TClonesArray *PHUrReader::fEvent = 0;
26
27PHUrReader::PHUrReader(TString inputfile) {
28
29 Input(inputfile);
30 collisionsIn .reserve(20000) ; // reserve enough space to avoid copying of objects during dynamic growing
31 collisionsOut.reserve(20000); // reserve enough space to avoid copying of objects during dynamic growing
32 collisions .reserve(2000); // reserve enough space to avoid copying of objects during dynamic growing
33 printEvtHeader = kTRUE;
34 printColHeader = kFALSE;
35 printParticle = kFALSE;
36 fEvtCt = 0;
37 fOutputAscii = NULL;
38 outputNonDiLeptons = kFALSE;
39 outputLeptons = kTRUE;
40 outputFreezeoutDiLeptons = kFALSE;
41 fNumMax = 10000;
42 /*
43 Taken from Unigen:
44 UrQMD PIDs are in fact composite - a particle is fully defined by the
45 type specifier (ityp), the charge (ichg) and in case of baryons, the
46 third component isospin (iso3; ignored by us at present). For
47 simplicity, our conversion tables collapse these into a single number
48 as follows:
49 - propagate the sign of ityp (particle-antiparticle distinction for
50 baryons, strangeness-antistrangeness distinction for mesons) to that
51 of the stored value;
52 - shift the ichg range from -2..2 (UrQMD does not support nuclear
53 fragments other than protons and neutrons so all particles it
54 produces fall in this range) to 0..4 to make sure it doesn't
55 interfere with the above;
56 - multiply shifted charge by +/-1000 and add it to type. The latter
57 is guaranteed to be smaller than 1000 (baryon types are one- or
58 two-digit, meson types three-digit) so that way no ambiguities
59 occur.
60 Int_t id;
61 if (ityp >= 0)
62 id = 1000 * (ichg + 2) + ityp;
63 else
64 id = -1000 * (ichg + 2) + ityp;
65 */
66 //
67 pdg_param = 0;
68 pid_param = 0;
69
70 mUrqmdToPdg[1017 ] = 1114;
71 mUrqmdToPdg[1018 ] =31114;
72 mUrqmdToPdg[1019 ] =1112 ;
73 mUrqmdToPdg[1020 ] =11114;
74 mUrqmdToPdg[1021 ] =11112;
75 mUrqmdToPdg[1022 ] =1116 ;
76 mUrqmdToPdg[1023 ] =21112;
77 mUrqmdToPdg[1024 ] =21114;
78 mUrqmdToPdg[1025 ] =11116;
79 mUrqmdToPdg[1026 ] =1118 ;
80 mUrqmdToPdg[1040 ] =3112 ;
81 mUrqmdToPdg[1041 ] =3114 ;
82 mUrqmdToPdg[1042 ] =13112;
83 mUrqmdToPdg[1043 ] =13114;
84 mUrqmdToPdg[1044 ] =23112;
85 mUrqmdToPdg[1045 ] =3116 ;
86 mUrqmdToPdg[1046 ] =13116;
87 mUrqmdToPdg[1047 ] =23114;
88 mUrqmdToPdg[1048 ] =3118 ;
89 mUrqmdToPdg[1049 ] =3312 ;
90 mUrqmdToPdg[1050 ] =3314 ;
91 mUrqmdToPdg[1051 ] =23314;
92 mUrqmdToPdg[1052 ] =13314;
93 mUrqmdToPdg[1053 ] =33314;
94 mUrqmdToPdg[1054 ] =13316;
95 mUrqmdToPdg[1055 ] =3334 ;
96 mUrqmdToPdg[1101 ] =-211 ;
97 mUrqmdToPdg[1104 ] =-213 ;
98 mUrqmdToPdg[1111 ] =-9000211;
99 mUrqmdToPdg[1114 ] =-20213 ;
100 mUrqmdToPdg[1118 ] =-215 ;
101 mUrqmdToPdg[1122 ] =-10213 ;
102 mUrqmdToPdg[1126 ] =-100213;
103 mUrqmdToPdg[1130 ] =-30213 ;
104 mUrqmdToPdg[2001 ] =2112 ;
105 mUrqmdToPdg[2002 ] =12112 ;
106 mUrqmdToPdg[2003 ] =1214 ;
107 mUrqmdToPdg[2004 ] =22112 ;
108 mUrqmdToPdg[2005 ] =32112 ;
109 mUrqmdToPdg[2006 ] =2116 ;
110 mUrqmdToPdg[2007 ] =12116 ;
111 mUrqmdToPdg[2008 ] =21214 ;
112 mUrqmdToPdg[2009 ] =42112 ;
113 mUrqmdToPdg[2010 ] =31214 ;
114 mUrqmdToPdg[2011 ] =41214 ;
115 mUrqmdToPdg[2012 ] =12118 ;
116 mUrqmdToPdg[2013 ] =52114 ;
117 mUrqmdToPdg[2016 ] =100012110;
118 mUrqmdToPdg[2017 ] =2114 ;
119 mUrqmdToPdg[2018 ] =32114 ;
120 mUrqmdToPdg[2019 ] =1212 ;
121 mUrqmdToPdg[2020 ] =12114 ;
122 mUrqmdToPdg[2021 ] =11212 ;
123 mUrqmdToPdg[2022 ] =1216 ;
124 mUrqmdToPdg[2023 ] =21212 ;
125 mUrqmdToPdg[2024 ] =22114 ;
126 mUrqmdToPdg[2025 ] =11216 ;
127 mUrqmdToPdg[2026 ] =2118 ;
128 mUrqmdToPdg[2027 ] =3122 ;
129 mUrqmdToPdg[2028 ] =13122 ;
130 mUrqmdToPdg[2029 ] =3124 ;
131 mUrqmdToPdg[2030 ] =23122 ;
132 mUrqmdToPdg[2031 ] =33122 ;
133 mUrqmdToPdg[2032 ] =13124 ;
134 mUrqmdToPdg[2033 ] =43122 ;
135 mUrqmdToPdg[2034 ] =53122 ;
136 mUrqmdToPdg[2035 ] =3126 ;
137 mUrqmdToPdg[2036 ] =13126 ;
138 mUrqmdToPdg[2037 ] =23124 ;
139 mUrqmdToPdg[2038 ] =3128 ;
140 mUrqmdToPdg[2039 ] =23126 ;
141 mUrqmdToPdg[2040 ] =3212 ;
142 mUrqmdToPdg[2041 ] =3214 ;
143 mUrqmdToPdg[2042 ] =13212 ;
144 mUrqmdToPdg[2043 ] =13214 ;
145 mUrqmdToPdg[2044 ] =23212 ;
146 mUrqmdToPdg[2045 ] =3216 ;
147 mUrqmdToPdg[2046 ] =13216 ;
148 mUrqmdToPdg[2047 ] =23214 ;
149 mUrqmdToPdg[2048 ] =3218 ;
150 mUrqmdToPdg[2049 ] =3322 ;
151 mUrqmdToPdg[2050 ] =3324 ;
152 mUrqmdToPdg[2051 ] =23324 ;
153 mUrqmdToPdg[2052 ] =13324 ;
154 mUrqmdToPdg[2053 ] =33324 ;
155 mUrqmdToPdg[2054 ] =13326 ;
156 mUrqmdToPdg[2100 ] =22 ;
157 mUrqmdToPdg[2101 ] =111 ;
158 mUrqmdToPdg[2102 ] =221 ;
159 mUrqmdToPdg[2103 ] =223 ;
160 mUrqmdToPdg[2104 ] =113 ;
161 mUrqmdToPdg[2105 ] =9000221 ;
162 mUrqmdToPdg[2106 ] =311 ;
163 mUrqmdToPdg[2107 ] =331 ;
164 mUrqmdToPdg[2108 ] =313 ;
165 mUrqmdToPdg[2109 ] =333 ;
166 mUrqmdToPdg[2110 ] =333 ;
167 mUrqmdToPdg[2111 ] =9000111 ;
168 mUrqmdToPdg[2112 ] =10221 ;
169 mUrqmdToPdg[2113 ] =20313 ;
170 mUrqmdToPdg[2114 ] =20113 ;
171 mUrqmdToPdg[2115 ] =20223 ;
172 mUrqmdToPdg[2116 ] =20333 ;
173 mUrqmdToPdg[2117 ] =315 ;
174 mUrqmdToPdg[2118 ] =115 ;
175 mUrqmdToPdg[2119 ] =225 ;
176 mUrqmdToPdg[2120 ] =335 ;
177 mUrqmdToPdg[2121 ] =10313 ;
178 mUrqmdToPdg[2122 ] =10113 ;
179 mUrqmdToPdg[2123 ] =10223 ;
180 mUrqmdToPdg[2124 ] =10333 ;
181 mUrqmdToPdg[2125 ] =100313 ;
182 mUrqmdToPdg[2126 ] =100113 ;
183 mUrqmdToPdg[2127 ] =100223 ;
184 mUrqmdToPdg[2128 ] =100333 ;
185 mUrqmdToPdg[2129 ] =30313 ;
186 mUrqmdToPdg[2130 ] =30113 ;
187 mUrqmdToPdg[2131 ] =30223 ;
188 mUrqmdToPdg[2132 ] =337 ;
189 mUrqmdToPdg[3001 ] =2212 ;
190 mUrqmdToPdg[3002 ] =12212 ;
191 mUrqmdToPdg[3003 ] =2124 ;
192 mUrqmdToPdg[3004 ] =22212 ;
193 mUrqmdToPdg[3005 ] =32212 ;
194 mUrqmdToPdg[3006 ] =2216 ;
195 mUrqmdToPdg[3007 ] =12216 ;
196 mUrqmdToPdg[3008 ] =22124 ;
197 mUrqmdToPdg[3009 ] =42212 ;
198 mUrqmdToPdg[3010 ] =32124 ;
199 mUrqmdToPdg[3011 ] =42124 ;
200 mUrqmdToPdg[3012 ] =12218 ;
201 mUrqmdToPdg[3013 ] =52214 ;
202 mUrqmdToPdg[3016 ] =100012210;
203 mUrqmdToPdg[3017 ] =2214 ;
204 mUrqmdToPdg[3018 ] =32214 ;
205 mUrqmdToPdg[3019 ] =2122 ;
206 mUrqmdToPdg[3020 ] =12214 ;
207 mUrqmdToPdg[3021 ] =12122 ;
208 mUrqmdToPdg[3022 ] =2126 ;
209 mUrqmdToPdg[3023 ] =22122 ;
210 mUrqmdToPdg[3024 ] =22214 ;
211 mUrqmdToPdg[3025 ] =12126 ;
212 mUrqmdToPdg[3026 ] =2218 ;
213 mUrqmdToPdg[3040 ] =3222 ;
214 mUrqmdToPdg[3041 ] =3224 ;
215 mUrqmdToPdg[3042 ] =13222 ;
216 mUrqmdToPdg[3043 ] =13224 ;
217 mUrqmdToPdg[3044 ] =23222 ;
218 mUrqmdToPdg[3045 ] =3226 ;
219 mUrqmdToPdg[3046 ] =13226 ;
220 mUrqmdToPdg[3047 ] =23224 ;
221 mUrqmdToPdg[3048 ] =3228 ;
222 mUrqmdToPdg[3101 ] =211 ;
223 mUrqmdToPdg[3104 ] =213 ;
224 mUrqmdToPdg[3106 ] =321 ;
225 mUrqmdToPdg[3108 ] =323 ;
226 mUrqmdToPdg[3110 ] =10321 ;
227 mUrqmdToPdg[3111 ] =9000211 ;
228 mUrqmdToPdg[3113 ] =20323 ;
229 mUrqmdToPdg[3114 ] =20213 ;
230 mUrqmdToPdg[3117 ] =325 ;
231 mUrqmdToPdg[3118 ] =215 ;
232 mUrqmdToPdg[3121 ] =10323 ;
233 mUrqmdToPdg[3122 ] =10213 ;
234 mUrqmdToPdg[3125 ] =100323 ;
235 mUrqmdToPdg[3126 ] =100213 ;
236 mUrqmdToPdg[3129 ] =30323 ;
237 mUrqmdToPdg[3130 ] =30213 ;
238 mUrqmdToPdg[4017 ] =2224 ;
239 mUrqmdToPdg[4018 ] =32224 ;
240 mUrqmdToPdg[4019 ] =2222 ;
241 mUrqmdToPdg[4020 ] =12224 ;
242 mUrqmdToPdg[4021 ] =12222 ;
243 mUrqmdToPdg[4022 ] =2226 ;
244 mUrqmdToPdg[4023 ] =22222 ;
245 mUrqmdToPdg[4024 ] =22224 ;
246 mUrqmdToPdg[4025 ] =12226 ;
247 mUrqmdToPdg[4026 ] =2228 ;
248 mUrqmdToPdg[-3055 ] = -3334 ;
249 mUrqmdToPdg[-3054 ] = -13316 ;
250 mUrqmdToPdg[-3053 ] = -33314 ;
251 mUrqmdToPdg[-3052 ] = -13314 ;
252 mUrqmdToPdg[-3051 ] = -23314 ;
253 mUrqmdToPdg[-3050 ] = -3314 ;
254 mUrqmdToPdg[-3049 ] = -3312 ;
255 mUrqmdToPdg[-3048 ] = -3118 ;
256 mUrqmdToPdg[-3047 ] = -23114 ;
257 mUrqmdToPdg[-3046 ] = -13116 ;
258 mUrqmdToPdg[-3045 ] = -3116 ;
259 mUrqmdToPdg[-3044 ] = -23112 ;
260 mUrqmdToPdg[-3043 ] = -13114 ;
261 mUrqmdToPdg[-3042 ] = -13112 ;
262 mUrqmdToPdg[-3041 ] = -3114 ;
263 mUrqmdToPdg[-3040 ] = -3112 ;
264 mUrqmdToPdg[-3026 ] = -1118 ;
265 mUrqmdToPdg[-3025 ] = -11116 ;
266 mUrqmdToPdg[-3024 ] = -21114 ;
267 mUrqmdToPdg[-3023 ] = -21112 ;
268 mUrqmdToPdg[-3022 ] = -1116 ;
269 mUrqmdToPdg[-3021 ] = -11112 ;
270 mUrqmdToPdg[-3020 ] = -11114 ;
271 mUrqmdToPdg[-3019 ] = -1112 ;
272 mUrqmdToPdg[-3018 ] = -31114 ;
273 mUrqmdToPdg[-3017 ] = -1114 ;
274 mUrqmdToPdg[-2129 ] = -30313 ;
275 mUrqmdToPdg[-2125 ] = -100313 ;
276 mUrqmdToPdg[-2121 ] = -10313 ;
277 mUrqmdToPdg[-2117 ] = -315 ;
278 mUrqmdToPdg[-2113 ] = -20313 ;
279 mUrqmdToPdg[-2110 ] = -10311 ;
280 mUrqmdToPdg[-2108 ] = -313 ;
281 mUrqmdToPdg[-2106 ] = -311 ;
282 mUrqmdToPdg[-2055 ] = -3334 ;
283 mUrqmdToPdg[-2054 ] = -13326 ;
284 mUrqmdToPdg[-2053 ] = -33324 ;
285 mUrqmdToPdg[-2052 ] = -13324 ;
286 mUrqmdToPdg[-2051 ] = -23324 ;
287 mUrqmdToPdg[-2050 ] = -3324 ;
288 mUrqmdToPdg[-2049 ] = -3322 ;
289 mUrqmdToPdg[-2048 ] = -3218 ;
290 mUrqmdToPdg[-2047 ] = -23214 ;
291 mUrqmdToPdg[-2046 ] = -13216 ;
292 mUrqmdToPdg[-2045 ] = -3216 ;
293 mUrqmdToPdg[-2044 ] = -23212 ;
294 mUrqmdToPdg[-2043 ] = -13214 ;
295 mUrqmdToPdg[-2042 ] = -13212 ;
296 mUrqmdToPdg[-2041 ] = -3214 ;
297 mUrqmdToPdg[-2040 ] = -3212 ;
298 mUrqmdToPdg[-2039 ] = -23126 ;
299 mUrqmdToPdg[-2038 ] = -3128 ;
300 mUrqmdToPdg[-2037 ] = -23124 ;
301 mUrqmdToPdg[-2036 ] = -13126 ;
302 mUrqmdToPdg[-2035 ] = -3126 ;
303 mUrqmdToPdg[-2034 ] = -53122 ;
304 mUrqmdToPdg[-2033 ] = -43122 ;
305 mUrqmdToPdg[-2032 ] = -13124 ;
306 mUrqmdToPdg[-2031 ] = -33122 ;
307 mUrqmdToPdg[-2030 ] = -23122 ;
308 mUrqmdToPdg[-2029 ] = -3124 ;
309 mUrqmdToPdg[-2028 ] = -13122 ;
310 mUrqmdToPdg[-2027 ] = -3122 ;
311 mUrqmdToPdg[-2026 ] = -2118 ;
312 mUrqmdToPdg[-2025 ] = -11216 ;
313 mUrqmdToPdg[-2024 ] = -22114 ;
314 mUrqmdToPdg[-2023 ] = -21212 ;
315 mUrqmdToPdg[-2022 ] = -1216 ;
316 mUrqmdToPdg[-2021 ] = -11212 ;
317 mUrqmdToPdg[-2020 ] = -12114 ;
318 mUrqmdToPdg[-2019 ] = -1212 ;
319 mUrqmdToPdg[-2018 ] = -32114 ;
320 mUrqmdToPdg[-2017 ] = -2114 ;
321 mUrqmdToPdg[-2016 ] = -100012110;
322 mUrqmdToPdg[-2013 ] = -52114 ;
323 mUrqmdToPdg[-2012 ] = -12118 ;
324 mUrqmdToPdg[-2011 ] = -41214 ;
325 mUrqmdToPdg[-2010 ] = -31214 ;
326 mUrqmdToPdg[-2009 ] = -42112 ;
327 mUrqmdToPdg[-2008 ] = -21214 ;
328 mUrqmdToPdg[-2007 ] = -12116 ;
329 mUrqmdToPdg[-2006 ] = -2116 ;
330 mUrqmdToPdg[-2005 ] = -32112 ;
331 mUrqmdToPdg[-2004 ] = -22112 ;
332 mUrqmdToPdg[-2003 ] = -1214 ;
333 mUrqmdToPdg[-2002 ] = -12112 ;
334 mUrqmdToPdg[-2001 ] = -2112 ;
335 mUrqmdToPdg[-1129 ] = -30323 ;
336 mUrqmdToPdg[-1125 ] = -100323 ;
337 mUrqmdToPdg[-1121 ] = -10323 ;
338 mUrqmdToPdg[-1117 ] = -325 ;
339 mUrqmdToPdg[-1113 ] = -20323 ;
340 mUrqmdToPdg[-1110 ] = -10321 ;
341 mUrqmdToPdg[-1108 ] = -323 ;
342 mUrqmdToPdg[-1106 ] = -321 ;
343 mUrqmdToPdg[-1048 ] = -3228 ;
344 mUrqmdToPdg[-1047 ] = -23224 ;
345 mUrqmdToPdg[-1046 ] = -13226 ;
346 mUrqmdToPdg[-1045 ] = -3226 ;
347 mUrqmdToPdg[-1044 ] = -23222 ;
348 mUrqmdToPdg[-1043 ] = -13224 ;
349 mUrqmdToPdg[-1042 ] = -13222 ;
350 mUrqmdToPdg[-1041 ] = -3224 ;
351 mUrqmdToPdg[-1040 ] = -3222 ;
352 mUrqmdToPdg[-1026 ] = -2218 ;
353 mUrqmdToPdg[-1025 ] = -12126 ;
354 mUrqmdToPdg[-1024 ] = -22214 ;
355 mUrqmdToPdg[-1023 ] = -22122 ;
356 mUrqmdToPdg[-1022 ] = -2126 ;
357 mUrqmdToPdg[-1021 ] = -12122 ;
358 mUrqmdToPdg[-1020 ] = -12214 ;
359 mUrqmdToPdg[-1019 ] = -2122 ;
360 mUrqmdToPdg[-1018 ] = -32214 ;
361 mUrqmdToPdg[-1017 ] = -2214 ;
362 mUrqmdToPdg[-1016 ] = -100012210;
363 mUrqmdToPdg[-1013 ] = -52214 ;
364 mUrqmdToPdg[-1012 ] = -12218 ;
365 mUrqmdToPdg[-1011 ] = -42124 ;
366 mUrqmdToPdg[-1010 ] = -32124 ;
367 mUrqmdToPdg[-1009 ] = -42212 ;
368 mUrqmdToPdg[-1008 ] = -22124 ;
369 mUrqmdToPdg[-1007 ] = -12216 ;
370 mUrqmdToPdg[-1006 ] = -2216 ;
371 mUrqmdToPdg[-1005 ] = -32212 ;
372 mUrqmdToPdg[-1004 ] = -22212 ;
373 mUrqmdToPdg[-1003 ] = -2124 ;
374 mUrqmdToPdg[-1002 ] = -12212 ;
375 mUrqmdToPdg[-1001 ] = -2212 ;
376 mUrqmdToPdg[-26 ] =-2228 ;
377 mUrqmdToPdg[-25 ] =-12226 ;
378 mUrqmdToPdg[-24 ] =-22224 ;
379 mUrqmdToPdg[-23 ] =-22222 ;
380 mUrqmdToPdg[-22 ] =-2226 ;
381 mUrqmdToPdg[-21 ] =-12222 ;
382 mUrqmdToPdg[-20 ] =-12224 ;
383 mUrqmdToPdg[-19 ] =-2222 ;
384 mUrqmdToPdg[-18 ] =-32224 ;
385 mUrqmdToPdg[-17 ] =-2224 ;
386
387
388 dilep.mUrqmdToPdg = &mUrqmdToPdg;
389
390 mUrqmdProcess[1] = "NN->ND";
391 mUrqmdProcess[2] = "NN->NN*";
392 mUrqmdProcess[3] = "NN->ND*";
393 mUrqmdProcess[4] = "NN->DD";
394 mUrqmdProcess[5] = "NN->DN*";
395 mUrqmdProcess[6] = "NN->DD*";
396 mUrqmdProcess[7] = "NN->N*N*,N*D*,D*D*";
397 mUrqmdProcess[8] = "ND->DD";
398 mUrqmdProcess[10] = "MB->B'";
399 mUrqmdProcess[11] = "MM->M'";
400 mUrqmdProcess[13] = "BB (but not pp,pn) elastic scattering";
401 mUrqmdProcess[14] = "inelastic scattering (no string excitation)";
402 mUrqmdProcess[15] = "BB->2 strings";
403 mUrqmdProcess[17] = "pn-elastic";
404 mUrqmdProcess[19] = "pp-elastic";
405 mUrqmdProcess[20] = "decay";
406 mUrqmdProcess[22] = "BBar elastic";
407 mUrqmdProcess[23] = "BBar annihilation->1 string";
408 mUrqmdProcess[24] = "BBar diffractive->2 strings";
409 mUrqmdProcess[26] = "MB elastic scattering";
410 mUrqmdProcess[27] = "MB,MM->1 string";
411 mUrqmdProcess[28] = "MB,MM->2 string";
412 mUrqmdProcess[30] = "ND->NN";
413 mUrqmdProcess[31] = "DD->DN";
414 mUrqmdProcess[32] = "DD->NN";
415 mUrqmdProcess[35] = "ND inelastic";
416 mUrqmdProcess[36] = "Danielewicz forward delay (MB->B')";
417 mUrqmdProcess[37] = "Danielewicz forward delay (MM->M')";
418 mUrqmdProcess[38] = "MM elastic scattering";
419 mUrqmdProcess[39] = "BBar inelastic scattering (no annihilation)";
420}
421
422PHUrReader::~PHUrReader() {
423 if(fInputAscii.is_open()) fInputAscii.close();
424}
425
426const Char_t *PHUrReader::GetUrQMDProcess(Int_t id) {
427 map<Int_t,TString>::iterator it = mUrqmdProcess.find(id);
428 if(it!=mUrqmdProcess.end()) return it->second.Data();
429 else Error("GetUrQMDProcess()", "Process %i unknown!", id);
430 return 0;
431}
432
433void PHUrReader::ClearVector(vector<vector<PHUrParticle> >& v) {
434
435 for(UInt_t i = 0; i < v.size(); i++){
436 v[i].clear();
437 }
438 v.clear();
439}
440
441void PHUrReader::Input(TString filename) {
442
443 if(filename.CompareTo("") != 0) {
444
445 if(gSystem->AccessPathName(filename.Data()) != 0){
446 Error("Input()", "File %s not does not exist!",filename.Data());
447 } else {
448 if(fInputAscii.is_open()) fInputAscii.close();
449 fInputName = filename;
450 fInputAscii.open(filename.Data());
451 }
452 } else Error("Input()", "File name not specified !");
453}
454
455void PHUrReader::Output(TString filename) {
456 if(filename.CompareTo("") != 0) {
457 fOutputName = filename;
458 } else Error("Output()", "File name not specified !");
459}
460
461void PHUrReader::SetMapID(Int_t Id_in, Int_t Id_out) {
462 // remap UrQMD id to another id
463 if(mapIDs.find(Id_in) == mapIDs.end()){
464 mapIDs[Id_in] = Id_out;
465 } else {
466 Error("SetMapID()", "ID=%i already mapped! Will be ignored ....", Id_in);
467 }
468}
469
470Bool_t PHUrReader::ReadEvent() {
471 // read function for UrQMD f15 files
472
473 if(!fInputAscii.is_open() || !fInputAscii.good()) return kFALSE;
474
475 if(fEvtCt == 0){
476 // only in first event. In next
477 // events the eventheader will be read
478 // by the previous call.
479 evtheader.Clear(NULL);
480 if(!evtheader.Read(fInputAscii)){
481 Error("ReadEvent()", "Could not read event header!");
482 return kFALSE;
483 }
484 } else {
485 evtheader = evtheaderCP;
486 }
487
488 if(evtheader.id == -1 || evtheader.id == 0) { // detected eventheader
489
490 if(printEvtHeader) { evtheader.Print(); }
491
492 ClearVector(collisionsIn); // clear collisions of previous event
493 ClearVector(collisionsOut);
494 collisions.clear();
495
496 while ( !fInputAscii.eof() ) {
497 colheader.Clear(NULL);
498 if(!colheader.Read(fInputAscii)) return kFALSE; // read next collision header
499
500 if(colheader.n_in == -1 || colheader.n_in == 0){ // detected eventheader
501 evtheaderCP.Copy(colheader); // move previous evtheader to current eventheader
502 fEvtCt++;
503 return kTRUE;
504 }
505
506 collisions.push_back(colheader);
507
508 if(printColHeader) colheader.Print();
509
510 //-------------------------------------------------------------------------------
511 // read in-going particles
512 vector< PHUrParticle > InParticles;
513
514 for(Int_t i = 0; i < colheader.n_in; i++) {
515 PHUrParticle part;
516 if(!part.Read(fInputAscii)) {
517 Error("ReadEvent()", "Could not read in-going particle!");
518 return kFALSE;
519 }
520
521 Int_t pdg = 1000;
522 Int_t id = part.pdg;
523 if(mUrqmdToPdg.find(id) != mUrqmdToPdg.end()) pdg = mUrqmdToPdg[id];
524 part.pdg = pdg;
525
526 if(printParticle) part.Print();
527 InParticles.push_back(part);
528 }
529 collisionsIn.push_back(InParticles);
530 //-------------------------------------------------------------------------------
531
532 //-------------------------------------------------------------------------------
533 // read out-going particles
534 vector< PHUrParticle > OutParticles;
535
536 for(Int_t i = 0; i < colheader.n_out; i++) {
537 PHUrParticle part;
538 if(!part.Read(fInputAscii)) {
539 Error("ReadEvent()", "Could not read out-going particle!");
540 return kFALSE;
541 }
542
543 Int_t pdg = 1000;
544 Int_t id = part.pdg;
545 if(mUrqmdToPdg.find(id) != mUrqmdToPdg.end()) pdg = mUrqmdToPdg[id];
546 part.pdg = pdg;
547
548 if(printParticle)part.Print();
549 OutParticles.push_back(part);
550 }
551 collisionsOut.push_back(OutParticles);
552 //-------------------------------------------------------------------------------
553
554 }
555 fEvtCt++; // last event only
556 return kTRUE;
557
558 } else {
559 Error("ReadEvent()", "This is not an event header!");
560 return kFALSE;
561 }
562
563}
564
565Bool_t PHUrReader::Modify(PParticle **mstack, int *, int *num, int stacksize) {
566
567 if(!ReadEvent()) return kFALSE;
568
569
570 //----------------------------------------------------------------------------------
571 // for the first call : create event structure
572 if(fEvtCt == 1) {
573
574 makeDistributionManager()->Exec("pdg:init");
575 pdg_param = makeDataBase()->GetParamInt("pdg");
576 if (pdg_param < 0) return kFALSE;
577 pid_param = makeDataBase()->GetParamInt("pid");
578
579 TTree *T = GetTree();
580
581 if(T) {
582 dilep.outputLeptons = outputLeptons;
583 dilep.reader = this;
584 if(outputLeptons) dilep.Output(fInputName,fOutputName);
585
586 if(stacksize == 0){
587 Error("Modify()", "FATAL ERROR : stacksize == 0!, Check your statcksize setiings for overrun!");
588 exit(1);
589 }
590
591 CreateParticleArray(T,stacksize );
592 CreateAddonArray (T,stacksize );
593
594 } else {
595 Error("Modify()", "retrieved NULL pointer for storage tree!");
596 return kFALSE;
597 }
598 }
599 //----------------------------------------------------------------------------------
600
601 //###########################################################################################
602 // loop over urqmd event
603 //c.. other
604 Int_t flag = -1 ;
605 Int_t flag_noLep = -1 ;
606 //c.. start analysis
607 //c.. loop over all interactions ww
608
609 UInt_t numberCol = collisions.size();
610 fNumPart = 0;
611 fNumPartAddon = 0;
612 fEvent->Clear();
613 fAdd ->Clear();
614
615 //----------------------------------------------------------------------------------
616 // map the occurence of one particle by index (for expample elastic scattering does not
617 // change the index of the particle, but momentum)
618
619
620 vector<PHUrParticle*> particles;
621 map <Int_t,PHUrParticle*> mLastParticle;
622
623 vector<Int_t> particlesCt;
624 particles .resize(20000,0);
625 particlesCt.resize(20000,0);
626
627 Int_t nparticle = 0;
628 for(UInt_t icol = 0; icol < numberCol; icol++) { // collision
629 for(UInt_t iin = 0; iin < collisionsIn[icol].size(); iin++) { // in-going particles
630 PHUrParticle *particle = &collisionsIn [icol][iin];
631 if(particles[particle->ind] == 0) {
632 particles[particle->ind] = particle;
633 nparticle++;
634 particle->first = 0;
635 mLastParticle[particle->ind] = particle;
636 } else {
637 mLastParticle[particle->ind] = particle;
638 }
639 particlesCt[particle->ind]++;
640 particle->instance = particlesCt[particle->ind];
641
642 particles[particle->ind]->listInCollIndex.push_back(icol);
643 }
644
645 for(UInt_t iout = 0; iout < collisionsOut[icol].size(); iout++) { // out-going particles
646 PHUrParticle *particle = &collisionsOut [icol][iout];
647 if(particles[particle->ind] == 0) {
648 particles[particle->ind] = particle;
649 nparticle++;
650 particlesCt[particle->ind]++;
651 particle->first = 1;
652 mLastParticle[particle->ind] = particle;
653 } else {
654 mLastParticle[particle->ind] = particle;
655 }
656 particle->instance = particlesCt[particle->ind];
657 particles[particle->ind]->listOutCollIndex.push_back(icol);
658 }
659 }
660
661 for(map< Int_t, PHUrParticle*>::iterator iter = mLastParticle.begin(); iter != mLastParticle.end(); ++iter ) {
662 if((*iter).second ) (*iter).second->last=1;
663 }
664
665 particles .erase(particles .begin()+nparticle, particles.end());
666 particlesCt.erase(particlesCt.begin()+nparticle,particlesCt.end());
667
668 //----------------------------------------------------------------------------------
669
670 //----------------------------------------------------------------------------------
671 if(outputNonDiLeptons == 2) {
672 Int_t ct = 0;
673 for(UInt_t icol = 0; icol < numberCol; icol++) { // collision
674 for(UInt_t iin = 0; iin < collisionsIn[icol].size(); iin++) { // in-going particles
675 PHUrParticle &particle = collisionsIn [icol][iin];
676 PHUrCollisionHeader &collision = collisions[icol];
677
678 if(particle.id == dilep.omega ||
679 particle.id == dilep.delta || // omega || delta || rho || phi
680 particle.id == dilep.rho ||
681 particle.id == dilep.phi ||
682 particle.id == dilep.pion ||
683 particle.id == dilep.etaprime ||
684 particle.id == dilep.eta
685 ) { // dilepton interest
686 continue;
687 }
688
689 if(ct == fNumMax-50) {
690 continue;
691 }
692
693 Int_t index = 0;
694 PParticle *pMoth = CreateParticle(index);
695 PHUrAddon *pMothAdd = CreateAddon(index);
696
697 pMoth->Reset(particle.pdg, particle.Px(), particle.Py(), particle.Pz(), particle.E(), 1.);
698
699 pMothAdd->t_cre = particle.t;
700 pMothAdd->t_abs = particle.t;
701 pMothAdd->dens_cre = collision.baryon_density;
702 pMothAdd->dens_abs = collision.baryon_density;
703 pMothAdd->fr_cre = particle.fr;
704 pMothAdd->fr_abs = particle.fr;
705 pMothAdd->stable = 0;
706 pMothAdd->n_in = collision.n_in;
707 pMothAdd->n_out = collision.n_out;
708 pMothAdd->process_id = collision.process_id;
709 pMothAdd->n_collision = collision.n_collision;
710 pMothAdd->instance = particle.instance;
711 pMothAdd->first = particle.first;
712 pMothAdd->last = particle.last;
713
714 }
715
716 for(UInt_t iout = 0; iout < collisionsOut[icol].size(); iout++) { // out-going particles
717
718 PHUrParticle &particle = collisionsOut [icol][iout];
719 PHUrCollisionHeader &collision = collisions[icol];
720
721 if(particle.id == dilep.omega ||
722 particle.id == dilep.delta || // omega || delta || rho || phi
723 particle.id == dilep.rho ||
724 particle.id == dilep.phi ||
725 particle.id == dilep.pion ||
726 particle.id == dilep.etaprime ||
727 particle.id == dilep.eta
728 ) { // dilepton interest
729 continue;
730 }
731
732 if(ct == fNumMax-50) {
733 continue;
734 }
735
736 Int_t index = 0;
737 PParticle *pMoth = CreateParticle(index);
738 PHUrAddon *pMothAdd = CreateAddon(index);
739
740 pMoth->Reset(particle.pdg, particle.Px(), particle.Py(), particle.Pz(), particle.E(), 1.);
741
742 pMothAdd->t_cre = particle.t;
743 pMothAdd->t_abs = particle .t;
744 pMothAdd->dens_cre = collision.baryon_density;
745 pMothAdd->dens_abs = collision.baryon_density;
746 pMothAdd->fr_cre = particle .fr;
747 pMothAdd->fr_abs = particle .fr;
748 pMothAdd->stable = 0;
749 pMothAdd->n_in = collision.n_in;
750 pMothAdd->n_out = collision.n_out;
751 pMothAdd->process_id = collision.process_id;
752 pMothAdd->n_collision = collision.n_collision;
753 pMothAdd->instance = particle.instance;
754 pMothAdd->first = particle.first;
755 pMothAdd->last = particle.last;
756
757 }
758 }
759 }
760 //----------------------------------------------------------------------------------
761
762
763
764 for(UInt_t icol = 0; icol < numberCol; icol++) { // collision
765 for(UInt_t iout = 0; iout < collisionsOut[icol].size(); iout++) { // out-going particles
766
767 //----------------------------------------------------------------------------------
768 // try to find the particle as incoming particle in the
769 // following collisions
770 for(UInt_t icol2 = icol+1; icol2 < numberCol; icol2++) { // following collisions
771 for(UInt_t iin = 0; iin < collisionsIn[icol2].size(); iin++) { // in-going particles
772
773 PHUrParticle &particleIn = collisionsIn [icol2][iin];
774 PHUrParticle &particleOut = collisionsOut[icol ][iout];
775
776 if (particleIn.id == particleOut.id && particleIn.Pz() == particleOut.Pz()) {
777 // same particle in following collision
778
779 // found the same particle
780 Bool_t dilepton = kFALSE;
781
782 if(particleIn.id == dilep.omega ||
783 particleIn.id == dilep.delta || // omega || delta || rho || phi
784 particleIn.id == dilep.rho ||
785 particleIn.id == dilep.phi) { // dilepton interest
786
787 dilep.evtheader = &evtheader;
788 dilep.particleIn = &particleIn;
789 dilep.particleOut = &particleOut;
790 dilep.collisionIn = &collisions[icol2];
791 dilep.collisionOut = &collisions[icol ];
792
793 if(!outputFreezeoutDiLeptons) {
794 dilep.Dilep();
795 } else
796 flag = 0; // not final state
797
798 dilepton = kTRUE;
799
800 } // if dilepton interest
801
802 if(particleIn.id == dilep.pion ||
803 particleIn.id == dilep.etaprime ||
804 particleIn.id == dilep.eta) { // pion || etaprime || eta
805 flag = 0;
806 dilepton = kTRUE;
807 }
808
809 if(!dilepton) flag_noLep = 0;
810
811 } // if same particle in following collision
812 } // for iin
813 } // for col2
814 //----------------------------------------------------------------------------------
815
816 PHUrParticle &particleOut = collisionsOut[icol ][iout];
817
818 //----------------------------------------------------------------------------------
819 // in case of freeze out dileptons we have to
820 // take care
821 if(particleOut.id == dilep.omega ||
822 particleOut.id == dilep.delta || // omega || delta || rho || phi
823 particleOut.id == dilep.rho ||
824 particleOut.id == dilep.phi) { // dilepton interest
825
826 if(outputFreezeoutDiLeptons && flag != 0) {
827
828 dilep.evtheader = &evtheader;
829 dilep.particleIn = &particleOut;
830 dilep.particleOut = &particleOut;
831 dilep.collisionIn = &collisions[icol];
832 dilep.collisionOut = &collisions[icol];
833
834 dilep.Dilep();
835 }
836 }
837 //----------------------------------------------------------------------------------
838
839
840 if(particleOut.id == dilep.pion ||
841 particleOut.id == dilep.etaprime ||
842 particleOut.id == dilep.eta) { // pion || etaprime || eta
843
844 if(flag != 0) { // particle has not undergone any following collsion (last state)
845
846 dilep.evtheader = &evtheader;
847 dilep.particleIn = &particleOut;
848 dilep.particleOut = &particleOut;
849 dilep.collisionIn = &collisions[icol];
850 dilep.collisionOut = &collisions[icol];
851
852 dilep.Dilep();
853 } // if flag != 0
854 } else {
855
856 if( !(particleOut.id == dilep.omega ||
857 particleOut.id == dilep.delta || // omega || delta || rho || phi
858 particleOut.id == dilep.rho ||
859 particleOut.id == dilep.phi )
860
861 ) {
862
863 if (flag_noLep != 0) { // final out non lepton
864
865 if (outputNonDiLeptons == 1) {
866
867 PHUrCollisionHeader &collisionOut = collisions[icol];
868 Int_t index = 0;
869 PParticle *pMoth = CreateParticle(index);
870 PHUrAddon *pMothAdd = CreateAddon(index);
871
872 pMoth->Reset(particleOut.pdg, particleOut.Px(), particleOut.Py(),
873 particleOut.Pz(), particleOut.E(), 1.);
874
875 pMothAdd->t_cre = particleOut.t;
876 pMothAdd->t_abs = particleOut .t;
877 pMothAdd->dens_cre = collisionOut.baryon_density;
878 pMothAdd->dens_abs = collisionOut.baryon_density;
879 pMothAdd->fr_cre = particleOut .fr;
880 pMothAdd->fr_abs = particleOut .fr;
881 pMothAdd->stable = 2;
882 pMothAdd->n_in = collisionOut.n_in;
883 pMothAdd->n_out = collisionOut.n_out;
884 pMothAdd->process_id = collisionOut.process_id;
885 pMothAdd->n_collision = collisionOut.n_collision;
886 pMothAdd->instance = particleOut.instance;
887 pMothAdd->first = particleOut.first;
888 pMothAdd->last = particleOut.last;
889
890 }
891 }
892 }
893 }
894
895 flag = 1;
896 flag_noLep = 1;
897 } // for iout
898 } // for icol
899
900 // end loop over urqmd event
901 //###########################################################################################
902
903
904 Int_t nParticle = fEvent->GetEntries();
905 Int_t startInd = *num;
906 Int_t cur = startInd;
907 Int_t *i_result = 0;
908 for (Int_t i = 0; i < nParticle; i++) {
909
910 if (i == stacksize) {
911 Warning("Modify", "Stack size too small, increase '_system_particle_stacksize'");
912 return kTRUE;
913 }
914
915 PParticle *part = (PParticle*) fEvent->At(i);
916 Int_t id = part->ID(); // already urqmd -> pdg code
917
918 if (!makeDataBase()->GetParamInt (pdg_param, id, pid_param, &i_result)) {
919
920 if(find(funknownIDs.begin(),funknownIDs.end(),id) == funknownIDs.end()){
921 funknownIDs.push_back(id);
922 Warning("Modify()",
923 "Particle with PDG code %i is not defined in Pluto, further particles will get the 'dummy' pid", id);
924 }
925
926 } else {
927
928 Int_t id = *i_result;
929 if(mapIDs.find(*i_result) != mapIDs.end()){ id = mapIDs[*i_result]; }
930 part->SetID(id);
931
932 }
933
934 *(mstack[cur]) = *part;
935
936 cur++;
937 }
938
939 *num += cur-startInd;
940
941 return kTRUE;
942}
943
944
945void PHUrReader::CreateAddonArray(TTree *T, Int_t stacksize) {
946 // create Addon TClonesArray and add it to the PLUTO tree
947
948 fAdd = new TClonesArray("PHUrAddon", stacksize);
949 T->Branch("Addon", &fAdd, 5000, 99);
950 fNumMax = stacksize;
951}
952
953PHUrAddon *PHUrReader::CreateAddon(Int_t &index) {
954 // create an Addon object at index
955 if (fNumPartAddon == fNumMax){
956 Warning("CreateAddon()", "Reached max stacksize, no addon created!");
957 return 0;
958 }
959
960 new((*fAdd)[fNumPartAddon++]) PHUrAddon();
961 index = fNumPartAddon - 1;
962 return (PHUrAddon*)((*fAdd)[fNumPartAddon - 1]);
963}
964
965void PHUrReader::CreateParticleArray(TTree *, Int_t stacksize) {
966 // create UrQMD PParticle local TClonesArray
967 fEvent = new TClonesArray( "PParticle", stacksize);
968 //T->Branch("Event", &fEvent,8000,99);
969 fNumMax = stacksize;
970}
971
972PParticle *PHUrReader::CreateParticle(Int_t &index) {
973 // create PParticle object in local event at index
974
975 if (fNumPart == fNumMax){
976 Warning("CreateAddon()", "Reached max stacksize, no addon created!");
977 return 0;
978 }
979
980 new((*fEvent)[fNumPart++]) PParticle();
981 index = fNumPart - 1;
982 return (PParticle*)((*fEvent)[fNumPart - 1]);
983}
984
985Int_t PHUrReader::UserAnalysis(PParticle**,Int_t npart) {
986
987 cout << "userAna nparticle = " << npart << endl;
988
989 return 1;
990}
991
T
Float_t T
Definition NiNi_2AGeV.C:7
ClassImp
ClassImp(PBeamLineSimulation) PBeamLineSimulation
Definition PBeamLineSimulation.cc:22
makeDataBase
PDataBase * makeDataBase()
makeDistributionManager
PDistributionManager * makeDistributionManager()
PBulkInterface::GetTree
TTree * GetTree(void)
Definition PBulkInterface.h:53
PDataBase::GetParamInt
Int_t GetParamInt(const char *paramname, Int_t length=-1)
PDistributionManager::Exec
Bool_t Exec(const char *command)
PHUrAddon::fr_abs
TVector3 fr_abs
Definition PHUrAddon.h:29
PHUrAddon::n_in
Int_t n_in
Definition PHUrAddon.h:20
PHUrAddon::process_id
Int_t process_id
Definition PHUrAddon.h:22
PHUrAddon::stable
Int_t stable
Definition PHUrAddon.h:30
PHUrAddon::instance
Int_t instance
Definition PHUrAddon.h:31
PHUrAddon::dens_abs
Double_t dens_abs
Definition PHUrAddon.h:27
PHUrAddon::t_abs
Double_t t_abs
Definition PHUrAddon.h:25
PHUrAddon::fr_cre
TVector3 fr_cre
Definition PHUrAddon.h:28
PHUrAddon::dens_cre
Double_t dens_cre
Definition PHUrAddon.h:26
PHUrAddon::last
Int_t last
Definition PHUrAddon.h:33
PHUrAddon::first
Int_t first
Definition PHUrAddon.h:32
PHUrAddon::t_cre
Double_t t_cre
Definition PHUrAddon.h:24
PHUrAddon::n_out
Int_t n_out
Definition PHUrAddon.h:21
PHUrAddon::n_collision
Int_t n_collision
Definition PHUrAddon.h:23
PHUrEventHeader::Read
Bool_t Read(ifstream &in)
Definition PHUrEventHeader.h:29
PHUrEventHeader::Clear
void Clear(Option_t *)
Definition PHUrEventHeader.h:67
PHUrParticle::Print
void Print(Option_t *)
Definition PHUrParticle.h:116
PHUrParticle::fr
TVector3 fr
Definition PHUrParticle.h:23
PHUrParticle::Read
Bool_t Read(ifstream &in)
Definition PHUrParticle.h:42
PHUrReader::Output
void Output(TString filename)
Definition PHUrReader.cc:455
PHUrReader::Input
void Input(TString filename)
Definition PHUrReader.cc:441
PHUrReader::CreateParticle
PParticle * CreateParticle(Int_t &index)
Definition PHUrReader.cc:972
PHUrReader::CreateAddonArray
void CreateAddonArray(TTree *T, Int_t stacksize)
Definition PHUrReader.cc:945
PHUrReader::Modify
virtual bool Modify(PParticle **array, int *decay_done, int *num, int maxnum)
Definition PHUrReader.cc:565
PHUrReader::PHUrReader
PHUrReader(TString inputfile="")
Definition PHUrReader.cc:27
PHUrReader::UserAnalysis
static Int_t UserAnalysis(PParticle **, Int_t)
Definition PHUrReader.cc:985
PHUrReader::SetMapID
void SetMapID(Int_t Id_in, Int_t Id_out)
Definition PHUrReader.cc:461
PHUrReader::CreateAddon
PHUrAddon * CreateAddon(Int_t &index)
Definition PHUrReader.cc:953
PHUrReader::CreateParticleArray
void CreateParticleArray(TTree *T, Int_t stacksize)
Definition PHUrReader.cc:965
PParticle::ID
Int_t ID() const
Definition PParticle.h:80
PParticle::SetID
void SetID(const Int_t id)
Definition PParticle.h:91
PParticle::Reset
void Reset(const Int_t id, const TLorentzVector &v, const Double_t w=1.)
Definition PParticle.h:125
dilepton
PParticle * dilepton
Definition embeddedParticles.C:16
filename
TString filename("eta_dalitz_a")
pdg
PPDGPlugin * pdg
Definition PDistributionManager.mod.cc:2