]>
Commit | Line | Data |
---|---|---|
2e967919 | 1 | /****************************************************************************** |
2 | * T H E R M I N A T O R * | |
3 | * THERMal heavy-IoN generATOR * | |
4 | * version 1.0 * | |
5 | * * | |
6 | * Authors of the model: Wojciech Broniowski, Wojciech.Broniowski@ifj.edu.pl, * | |
7 | * Wojciech Florkowski, Wojciech.Florkowski@ifj.edu.pl * | |
8 | * Authors of the code: Adam Kisiel, kisiel@if.pw.edu.pl * | |
9 | * Tomasz Taluc, ttaluc@if.pw.edu.pl * | |
10 | * Code designers: Adam Kisiel, Tomasz Taluc, Wojciech Broniowski, * | |
11 | * Wojciech Florkowski * | |
12 | * * | |
13 | * For the detailed description of the program and furhter references * | |
14 | * to the description of the model plesase refer to: nucl-th/0504047, * | |
15 | * accessibile at: http://www.arxiv.org/nucl-th/0504047 * | |
16 | * * | |
17 | * Homepage: http://hirg.if.pw.edu.pl/en/therminator/ * | |
18 | * * | |
19 | * This code can be freely used and redistributed. However if you decide to * | |
20 | * make modifications to the code, please contact the authors, especially * | |
21 | * if you plan to publish the results obtained with such modified code. * | |
22 | * Any publication of results obtained using this code must include the * | |
23 | * reference to nucl-th/0504047 and the published version of it, when * | |
24 | * available. * | |
25 | * * | |
26 | *****************************************************************************/ | |
27 | #include <iostream> | |
28 | #include <TFile.h> | |
29 | #include <TH1D.h> | |
30 | #include "Particle.h" | |
31 | #include "Event.h" | |
32 | #include "Integrator.h" | |
33 | #include "ParticleDecayer.h" | |
34 | #include "DecayTable.h" | |
35 | #include "ReadPar.h" | |
36 | ||
37 | extern ReadPar* sRPInstance; | |
38 | ||
39 | Event::Event(ParticleDB *aDB, Integrator* aInteg) : | |
06b10510 | 40 | mCurIter(), |
2e967919 | 41 | mScanStarted(0) |
42 | { | |
43 | mDB = aDB; | |
44 | mInteg = aInteg; | |
45 | mRandom = new TRandom2(); | |
46 | #ifdef _ROOT_4_ | |
47 | mRandom->SetSeed2(31851, 14327); | |
48 | #else | |
49 | mRandom->SetSeed(31851); | |
50 | #endif | |
51 | mAverageMultiplicities.resize(mDB->GetParticleTypeCount()); | |
52 | mMultiplicities.resize(mDB->GetParticleTypeCount()); | |
53 | ReadParameters(); | |
54 | ReadMultiplicities(); | |
55 | ||
56 | ||
57 | } | |
58 | ||
59 | Event::~Event() | |
60 | { | |
61 | ||
62 | mParticles.clear(); | |
63 | delete mRandom; | |
64 | } | |
65 | ||
66 | void | |
67 | Event::GenerateEvent(int aSeed) | |
68 | { | |
69 | Particle **tPartBuf; | |
70 | ||
71 | mParticles.clear(); | |
72 | #ifdef _ROOT_4_ | |
73 | if (aSeed) mRandom->SetSeed2(aSeed, (aSeed*2) % (7*11*23*31)); | |
74 | #else | |
75 | if (aSeed) mRandom->SetSeed(aSeed); | |
76 | #endif | |
77 | GenerateMultiplicities(); | |
78 | for (unsigned int tIter=0; tIter<mAverageMultiplicities.size(); tIter++) { | |
79 | mInteg->Generate(mDB->GetParticleType(tIter), mMultiplicities[tIter], &tPartBuf); | |
80 | for (int tBufIter=0; tBufIter<mMultiplicities[tIter]; tBufIter++) { | |
81 | AddParticle(tPartBuf[tBufIter]); | |
82 | delete tPartBuf[tBufIter]; | |
83 | } | |
84 | free(tPartBuf); | |
85 | } | |
86 | } | |
87 | ||
88 | int | |
89 | Event::GetParticleCount() | |
90 | { | |
91 | return mParticles.size(); | |
92 | } | |
93 | ||
94 | Particle* | |
95 | Event::GetNextParticle() | |
96 | { | |
97 | if (mScanStarted) | |
98 | InitParticleScan(); | |
99 | if (mCurIter == mParticles.end()) | |
100 | return 0; | |
101 | ||
102 | Particle *tBuf; | |
103 | tBuf = &(*mCurIter); | |
104 | mCurIter++; | |
105 | ||
106 | return tBuf; | |
107 | } | |
108 | ||
109 | void | |
110 | Event::InitParticleScan() | |
111 | { | |
112 | mCurIter = mParticles.begin(); | |
113 | mScanStarted = 1; | |
114 | } | |
115 | ||
116 | Particle* | |
117 | Event::GetParticleOfCount(int aCount) | |
118 | { | |
119 | if (aCount == 0) { | |
120 | InitParticleScan(); | |
121 | } | |
122 | else if (aCount == 1) { | |
123 | InitParticleScan(); | |
124 | mCurIter++; | |
125 | } | |
126 | else { | |
127 | mCurIter--; | |
128 | mCurIter++; | |
129 | mCurIter++; | |
130 | } | |
131 | ||
132 | if (mCurIter == mParticles.end()) | |
133 | return 0; | |
134 | ||
135 | Particle *tBuf; | |
136 | tBuf = &(*mCurIter); | |
137 | return tBuf; | |
138 | } | |
139 | ||
140 | void | |
141 | Event::WriteEvent(int nEvent) | |
142 | { | |
143 | ParticleListIter tPLIter; | |
144 | tPLIter = mParticles.begin(); | |
145 | int tCount = 0; | |
146 | ofstream os; | |
147 | ||
148 | if (nEvent == 0) { | |
149 | try { | |
150 | os.open(sRPInstance->getPar("EventOutputFile").Data()); | |
151 | } | |
152 | catch (STR e) { | |
153 | PRINT_MESSAGE("Very strange: No event output filename ???"); | |
154 | exit(0); | |
155 | } | |
156 | } | |
157 | else { | |
158 | try { | |
159 | os.open(sRPInstance->getPar("EventOutputFile").Data(), ios::app); | |
160 | } | |
161 | catch (STR e) { | |
162 | PRINT_MESSAGE("Very strange: No event output filename ???"); | |
163 | exit(0); | |
164 | } | |
165 | } | |
166 | ||
167 | ||
168 | ||
169 | if (nEvent == 0) | |
170 | { | |
171 | os << "Therminator text output\nall_particles_p_x\nTherminator " << endl; | |
172 | } | |
173 | os << nEvent+1 << " " << GetParticleCount() << " " << "0.000" << " " << "0.000" << endl; | |
174 | ||
175 | while (tPLIter != mParticles.end()) | |
176 | { | |
177 | os.flags(ios::right); | |
178 | os.precision(6); | |
179 | os.width(6); | |
180 | os << tCount << " "; | |
181 | (*tPLIter).WriteParticle(&os); | |
182 | os << endl; | |
183 | ||
184 | tCount++; | |
185 | tPLIter++; | |
186 | } | |
187 | } | |
188 | ||
189 | ||
190 | void | |
191 | Event::AddParticle(Particle* aParticle) | |
192 | { | |
193 | if (0) { | |
194 | ParticleListIter tPLIter; | |
195 | tPLIter = mParticles.begin(); | |
196 | while (tPLIter != mParticles.end()) | |
197 | { | |
198 | if (((*tPLIter).GetMass() == aParticle->GetMass()) && | |
199 | ((*tPLIter).GetI3() == aParticle->GetI3()) && | |
200 | ((*tPLIter).GetBarionN() == aParticle->GetBarionN()) && | |
d02d3abd | 201 | ((*tPLIter).GetStrangeness() == aParticle->GetStrangeness())) { |
202 | break; } | |
203 | if ((*tPLIter).GetMass() < aParticle->GetMass()) { | |
204 | break; } | |
205 | else if ((*tPLIter).GetMass() == aParticle->GetMass()) { | |
206 | if ((*tPLIter).GetI3() < aParticle->GetI3()) { | |
207 | break; } | |
208 | else if ((*tPLIter).GetI3() == aParticle->GetI3()) { | |
209 | if ((*tPLIter).GetBarionN() < aParticle->GetBarionN()) { | |
210 | break; } | |
211 | else if ((*tPLIter).GetBarionN() == aParticle->GetBarionN()) { | |
212 | if ((*tPLIter).GetStrangeness() < aParticle->GetStrangeness()) { | |
213 | break; } | |
214 | } | |
215 | } | |
216 | } | |
2e967919 | 217 | tPLIter++; |
218 | } | |
219 | mParticles.insert(tPLIter, *aParticle); | |
220 | } | |
221 | else | |
222 | mParticles.push_back(*aParticle); | |
223 | } | |
224 | ||
225 | void | |
226 | Event::ReadMultiplicities() | |
227 | { | |
228 | char tName[200]; | |
229 | double tMult; | |
230 | ||
231 | char *tHash; | |
13e7fe11 | 232 | int multsize; |
233 | char *tMultName = 0; | |
2e967919 | 234 | ifstream *fin = NULL; |
235 | ||
236 | tHash = mInteg->ParameterHash(); | |
13e7fe11 | 237 | multsize = strlen(mFOHSlocation.Data()) + 25 + strlen(tHash); |
238 | tMultName = (char *) malloc(sizeof(char) * multsize); | |
239 | ||
240 | if (!tMultName) { | |
241 | printf("Cannot allocate memory!\n"); | |
242 | exit(0); | |
243 | } | |
244 | ||
ae89e57a | 245 | if (mFOHSlocation != "") { |
246 | strcpy(tMultName, mFOHSlocation.Data()); | |
247 | strcat(tMultName, "/"); | |
248 | strcat(tMultName, "fmultiplicity_"); | |
249 | strcat(tMultName, tHash); | |
250 | strcat(tMultName, ".txt"); | |
251 | fin = new ifstream(tMultName); | |
252 | } | |
253 | else if (!((fin) && (fin->is_open()))) { | |
254 | strcpy(tMultName, "fmultiplicity_"); | |
255 | strcat(tMultName, tHash); | |
256 | strcat(tMultName, ".txt"); | |
257 | fin = new ifstream(tMultName); | |
258 | } | |
259 | ||
260 | // strcpy(tMultName, "fmultiplicity_"); | |
261 | // strcat(tMultName, tHash); | |
262 | // strcat(tMultName, ".txt"); | |
2e967919 | 263 | |
13e7fe11 | 264 | // fin = new ifstream(tMultName); |
2e967919 | 265 | if ((fin) && (fin->is_open())) { |
266 | PRINT_MESSAGE("Reading Multiplicities from " << tMultName); | |
267 | ||
268 | while (!fin->eof()) | |
269 | { | |
270 | (*fin) >> tName >> tMult; | |
271 | PRINT_DEBUG_2(tName << " " << mDB->GetParticleTypeIndex(strdup(tName)) << " " << tMult); | |
272 | mAverageMultiplicities[mDB->GetParticleTypeIndex(strdup(tName))] = tMult; | |
273 | } | |
274 | fin->close(); | |
275 | } | |
13e7fe11 | 276 | |
277 | delete tMultName; | |
2e967919 | 278 | } |
279 | ||
280 | void | |
281 | Event::GenerateMultiplicities() | |
282 | { | |
283 | for (unsigned int tIter=0; tIter<mAverageMultiplicities.size(); tIter++) | |
284 | mMultiplicities[tIter] = mRandom->Poisson(mAverageMultiplicities[tIter]); | |
285 | } | |
286 | ||
287 | void | |
288 | Event::DecayParticles() | |
289 | { | |
290 | Particle *tPart1, *tPart2, *tPart3, *tFather; | |
291 | ParticleDecayer* tDecayer; | |
292 | int tCount = 0; | |
293 | ||
294 | tDecayer = new ParticleDecayer(mDB); | |
295 | tDecayer->SeedSet(mRandom->Integer(100000000)); | |
296 | ||
297 | // InitParticleScan(); | |
298 | while ((tFather = GetParticleOfCount(tCount))) | |
299 | { | |
d02d3abd | 300 | if (tFather->GetParticleType()->GetGamma() >= 0.0) { |
2e967919 | 301 | if ((tFather->GetParticleType()->GetTable()) && (((DecayTable *) tFather->GetParticleType()->GetTable())->GetChannelCount()+1 > 0)) |
302 | { | |
303 | tDecayer->DecayParticle(tFather, &tPart1, &tPart2, &tPart3); | |
304 | #ifndef _RESCALE_CHANNELS_ | |
305 | if (!tPart1) { | |
306 | tCount++; | |
307 | continue; | |
308 | } | |
309 | #endif | |
310 | #ifdef _NO_THREE_BODY_DECAYS_ | |
311 | if (tPart3){ | |
312 | tCount++; | |
313 | continue; | |
314 | } | |
315 | #endif | |
316 | #ifdef _OMIT_TWO_BODY_ | |
317 | if (tPart1 && tPart2 && !tPart3){ | |
318 | tCount++; | |
319 | continue; | |
320 | } | |
321 | #endif | |
322 | ||
323 | tPart1->SetFather(tCount); | |
324 | tPart2->SetFather(tCount); | |
325 | if (tPart3) | |
326 | tPart3->SetFather(tCount); | |
327 | ||
328 | AddParticle(tPart1); | |
329 | AddParticle(tPart2); | |
330 | ||
331 | delete tPart1; | |
332 | delete tPart2; | |
333 | ||
334 | if (tPart3) { | |
335 | AddParticle(tPart3); | |
336 | delete tPart3; | |
337 | } | |
338 | } | |
339 | else | |
340 | { | |
341 | } | |
d02d3abd | 342 | } |
2e967919 | 343 | tCount++; |
344 | } | |
345 | ||
346 | delete tDecayer; | |
347 | } | |
348 | ||
349 | void | |
350 | Event::Randomize() | |
351 | { | |
352 | TDatime dat; | |
353 | ||
354 | #ifdef _ROOT_4_ | |
355 | mRandom->SetSeed2(dat.Get() / 2 * 3, dat.Get() / 11 * 9); | |
356 | #else | |
357 | mRandom->SetSeed(dat.Get() / 2 * 3); | |
358 | #endif | |
359 | } | |
360 | ||
361 | void | |
362 | Event::ReadParameters() | |
363 | { | |
364 | STR tEvFile; | |
365 | STR tUseNegativeBinomial; | |
366 | ||
367 | try { | |
368 | tEvFile = sRPInstance->getPar("EventOutputFile"); | |
369 | } | |
370 | catch (STR e) { | |
371 | PRINT_DEBUG_1("Event::ReadParameters - Caught exception " << e); | |
372 | PRINT_MESSAGE("Did not find one of the neccessary parameters in the parameters file."); | |
373 | exit(0); | |
374 | } | |
375 | try { | |
376 | tUseNegativeBinomial = sRPInstance->getPar("MultiplicityDistribution"); | |
377 | if (tUseNegativeBinomial.Contains("NegativeBinomial")) | |
378 | mNegBin = 1; | |
379 | else | |
380 | mNegBin = 0; | |
381 | } | |
382 | catch (STR e) { | |
383 | PRINT_MESSAGE("Using default multiplicty distribution: Poissonian"); | |
384 | mNegBin = 0; | |
385 | } | |
ae89e57a | 386 | try { |
387 | mFOHSlocation = sRPInstance->getPar("FOHSLocation"); | |
388 | } | |
389 | catch (STR e) { | |
390 | } | |
2e967919 | 391 | |
392 | } | |
393 |