Update master to aliroot
[u/mrichter/AliRoot.git] / TTherminator / Therminator / Event.cxx
CommitLineData
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
37extern ReadPar* sRPInstance;
38
39Event::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
59Event::~Event()
60{
61
62 mParticles.clear();
63 delete mRandom;
64}
65
66void
67Event::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
88int
89Event::GetParticleCount()
90{
91 return mParticles.size();
92}
93
94Particle*
95Event::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
109void
110Event::InitParticleScan()
111{
112 mCurIter = mParticles.begin();
113 mScanStarted = 1;
114}
115
116Particle*
117Event::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
140void
141Event::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
190void
191Event::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
225void
226Event::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;
09b7c66c 271 PRINT_DEBUG_2(tName << " " << mDB->GetParticleTypeIndex(tName) << " " << tMult);
d58024dc 272 mAverageMultiplicities[mDB->GetParticleTypeIndex(tName)] = tMult;
2e967919 273 }
274 fin->close();
275 }
13e7fe11 276
09b7c66c 277 delete tHash;
13e7fe11 278 delete tMultName;
2e967919 279}
280
281void
282Event::GenerateMultiplicities()
283{
284 for (unsigned int tIter=0; tIter<mAverageMultiplicities.size(); tIter++)
285 mMultiplicities[tIter] = mRandom->Poisson(mAverageMultiplicities[tIter]);
286}
287
288void
289Event::DecayParticles()
290{
f85eb564 291 Particle *tPart1=NULL, *tPart2=NULL, *tPart3=NULL, *tFather=NULL;
2e967919 292 ParticleDecayer* tDecayer;
293 int tCount = 0;
294
295 tDecayer = new ParticleDecayer(mDB);
296 tDecayer->SeedSet(mRandom->Integer(100000000));
297
298 // InitParticleScan();
299 while ((tFather = GetParticleOfCount(tCount)))
300 {
d02d3abd 301 if (tFather->GetParticleType()->GetGamma() >= 0.0) {
2e967919 302 if ((tFather->GetParticleType()->GetTable()) && (((DecayTable *) tFather->GetParticleType()->GetTable())->GetChannelCount()+1 > 0))
303 {
304 tDecayer->DecayParticle(tFather, &tPart1, &tPart2, &tPart3);
305#ifndef _RESCALE_CHANNELS_
a369fb20 306 if ((!tPart1) || (!tPart2)) {
2e967919 307 tCount++;
a369fb20 308 if (tPart1) delete tPart1;
309 if (tPart2) delete tPart2;
2e967919 310 continue;
311 }
312#endif
313#ifdef _NO_THREE_BODY_DECAYS_
314 if (tPart3){
315 tCount++;
6750cf94 316 if (tPart1) delete tPart1;
317 if (tPart2) delete tPart2;
2e967919 318 continue;
319 }
320#endif
321#ifdef _OMIT_TWO_BODY_
322 if (tPart1 && tPart2 && !tPart3){
323 tCount++;
6750cf94 324 delete tPart1;
325 delete tPart2;
2e967919 326 continue;
327 }
328#endif
329
330 tPart1->SetFather(tCount);
331 tPart2->SetFather(tCount);
332 if (tPart3)
333 tPart3->SetFather(tCount);
334
335 AddParticle(tPart1);
336 AddParticle(tPart2);
337
338 delete tPart1;
339 delete tPart2;
340
341 if (tPart3) {
342 AddParticle(tPart3);
343 delete tPart3;
344 }
345 }
346 else
347 {
348 }
d02d3abd 349 }
2e967919 350 tCount++;
351 }
352
353 delete tDecayer;
354}
355
356void
357Event::Randomize()
358{
359 TDatime dat;
360
361#ifdef _ROOT_4_
362 mRandom->SetSeed2(dat.Get() / 2 * 3, dat.Get() / 11 * 9);
363#else
364 mRandom->SetSeed(dat.Get() / 2 * 3);
365#endif
366}
367
368void
369Event::ReadParameters()
370{
371 STR tEvFile;
372 STR tUseNegativeBinomial;
373
374 try {
375 tEvFile = sRPInstance->getPar("EventOutputFile");
376 }
377 catch (STR e) {
378 PRINT_DEBUG_1("Event::ReadParameters - Caught exception " << e);
379 PRINT_MESSAGE("Did not find one of the neccessary parameters in the parameters file.");
380 exit(0);
381 }
382 try {
383 tUseNegativeBinomial = sRPInstance->getPar("MultiplicityDistribution");
384 if (tUseNegativeBinomial.Contains("NegativeBinomial"))
385 mNegBin = 1;
386 else
387 mNegBin = 0;
388 }
389 catch (STR e) {
390 PRINT_MESSAGE("Using default multiplicty distribution: Poissonian");
391 mNegBin = 0;
392 }
ae89e57a 393 try {
394 mFOHSlocation = sRPInstance->getPar("FOHSLocation");
395 }
396 catch (STR e) {
397 }
2e967919 398
399}
400