]>
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()) && | |
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 | tPLIter++; | |
215 | } | |
216 | mParticles.insert(tPLIter, *aParticle); | |
217 | } | |
218 | else | |
219 | mParticles.push_back(*aParticle); | |
220 | } | |
221 | ||
222 | void | |
223 | Event::ReadMultiplicities() | |
224 | { | |
225 | char tName[200]; | |
226 | double tMult; | |
227 | ||
228 | char *tHash; | |
13e7fe11 | 229 | int multsize; |
230 | char *tMultName = 0; | |
2e967919 | 231 | ifstream *fin = NULL; |
232 | ||
233 | tHash = mInteg->ParameterHash(); | |
13e7fe11 | 234 | multsize = strlen(mFOHSlocation.Data()) + 25 + strlen(tHash); |
235 | tMultName = (char *) malloc(sizeof(char) * multsize); | |
236 | ||
237 | if (!tMultName) { | |
238 | printf("Cannot allocate memory!\n"); | |
239 | exit(0); | |
240 | } | |
241 | ||
ae89e57a | 242 | if (mFOHSlocation != "") { |
243 | strcpy(tMultName, mFOHSlocation.Data()); | |
244 | strcat(tMultName, "/"); | |
245 | strcat(tMultName, "fmultiplicity_"); | |
246 | strcat(tMultName, tHash); | |
247 | strcat(tMultName, ".txt"); | |
248 | fin = new ifstream(tMultName); | |
249 | } | |
250 | else if (!((fin) && (fin->is_open()))) { | |
251 | strcpy(tMultName, "fmultiplicity_"); | |
252 | strcat(tMultName, tHash); | |
253 | strcat(tMultName, ".txt"); | |
254 | fin = new ifstream(tMultName); | |
255 | } | |
256 | ||
257 | // strcpy(tMultName, "fmultiplicity_"); | |
258 | // strcat(tMultName, tHash); | |
259 | // strcat(tMultName, ".txt"); | |
2e967919 | 260 | |
13e7fe11 | 261 | // fin = new ifstream(tMultName); |
2e967919 | 262 | if ((fin) && (fin->is_open())) { |
263 | PRINT_MESSAGE("Reading Multiplicities from " << tMultName); | |
264 | ||
265 | while (!fin->eof()) | |
266 | { | |
267 | (*fin) >> tName >> tMult; | |
268 | PRINT_DEBUG_2(tName << " " << mDB->GetParticleTypeIndex(strdup(tName)) << " " << tMult); | |
269 | mAverageMultiplicities[mDB->GetParticleTypeIndex(strdup(tName))] = tMult; | |
270 | } | |
271 | fin->close(); | |
272 | } | |
13e7fe11 | 273 | |
274 | delete tMultName; | |
2e967919 | 275 | } |
276 | ||
277 | void | |
278 | Event::GenerateMultiplicities() | |
279 | { | |
280 | for (unsigned int tIter=0; tIter<mAverageMultiplicities.size(); tIter++) | |
281 | mMultiplicities[tIter] = mRandom->Poisson(mAverageMultiplicities[tIter]); | |
282 | } | |
283 | ||
284 | void | |
285 | Event::DecayParticles() | |
286 | { | |
287 | Particle *tPart1, *tPart2, *tPart3, *tFather; | |
288 | ParticleDecayer* tDecayer; | |
289 | int tCount = 0; | |
290 | ||
291 | tDecayer = new ParticleDecayer(mDB); | |
292 | tDecayer->SeedSet(mRandom->Integer(100000000)); | |
293 | ||
294 | // InitParticleScan(); | |
295 | while ((tFather = GetParticleOfCount(tCount))) | |
296 | { | |
297 | if (tFather->GetParticleType()->GetGamma() >= 0.0) | |
298 | if ((tFather->GetParticleType()->GetTable()) && (((DecayTable *) tFather->GetParticleType()->GetTable())->GetChannelCount()+1 > 0)) | |
299 | { | |
300 | tDecayer->DecayParticle(tFather, &tPart1, &tPart2, &tPart3); | |
301 | #ifndef _RESCALE_CHANNELS_ | |
302 | if (!tPart1) { | |
303 | tCount++; | |
304 | continue; | |
305 | } | |
306 | #endif | |
307 | #ifdef _NO_THREE_BODY_DECAYS_ | |
308 | if (tPart3){ | |
309 | tCount++; | |
310 | continue; | |
311 | } | |
312 | #endif | |
313 | #ifdef _OMIT_TWO_BODY_ | |
314 | if (tPart1 && tPart2 && !tPart3){ | |
315 | tCount++; | |
316 | continue; | |
317 | } | |
318 | #endif | |
319 | ||
320 | tPart1->SetFather(tCount); | |
321 | tPart2->SetFather(tCount); | |
322 | if (tPart3) | |
323 | tPart3->SetFather(tCount); | |
324 | ||
325 | AddParticle(tPart1); | |
326 | AddParticle(tPart2); | |
327 | ||
328 | delete tPart1; | |
329 | delete tPart2; | |
330 | ||
331 | if (tPart3) { | |
332 | AddParticle(tPart3); | |
333 | delete tPart3; | |
334 | } | |
335 | } | |
336 | else | |
337 | { | |
338 | } | |
339 | tCount++; | |
340 | } | |
341 | ||
342 | delete tDecayer; | |
343 | } | |
344 | ||
345 | void | |
346 | Event::Randomize() | |
347 | { | |
348 | TDatime dat; | |
349 | ||
350 | #ifdef _ROOT_4_ | |
351 | mRandom->SetSeed2(dat.Get() / 2 * 3, dat.Get() / 11 * 9); | |
352 | #else | |
353 | mRandom->SetSeed(dat.Get() / 2 * 3); | |
354 | #endif | |
355 | } | |
356 | ||
357 | void | |
358 | Event::ReadParameters() | |
359 | { | |
360 | STR tEvFile; | |
361 | STR tUseNegativeBinomial; | |
362 | ||
363 | try { | |
364 | tEvFile = sRPInstance->getPar("EventOutputFile"); | |
365 | } | |
366 | catch (STR e) { | |
367 | PRINT_DEBUG_1("Event::ReadParameters - Caught exception " << e); | |
368 | PRINT_MESSAGE("Did not find one of the neccessary parameters in the parameters file."); | |
369 | exit(0); | |
370 | } | |
371 | try { | |
372 | tUseNegativeBinomial = sRPInstance->getPar("MultiplicityDistribution"); | |
373 | if (tUseNegativeBinomial.Contains("NegativeBinomial")) | |
374 | mNegBin = 1; | |
375 | else | |
376 | mNegBin = 0; | |
377 | } | |
378 | catch (STR e) { | |
379 | PRINT_MESSAGE("Using default multiplicty distribution: Poissonian"); | |
380 | mNegBin = 0; | |
381 | } | |
ae89e57a | 382 | try { |
383 | mFOHSlocation = sRPInstance->getPar("FOHSLocation"); | |
384 | } | |
385 | catch (STR e) { | |
386 | } | |
2e967919 | 387 | |
388 | } | |
389 |