Another attempt at Coverity report
[u/mrichter/AliRoot.git] / TTherminator / Therminator / Event.cxx
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) :
40   mCurIter(),
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             }
215           }
216         }
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;
232   int   multsize;
233   char *tMultName = 0;
234   ifstream *fin = NULL;
235   
236   tHash = mInteg->ParameterHash();
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
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");
263
264 //  fin = new ifstream(tMultName);
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(tName) << " " << tMult);
272       mAverageMultiplicities[mDB->GetParticleTypeIndex(tName)] = tMult;
273     }
274     fin->close();
275   }
276
277   delete tHash;
278   delete tMultName;
279 }
280
281 void           
282 Event::GenerateMultiplicities()
283 {
284   for (unsigned int tIter=0; tIter<mAverageMultiplicities.size(); tIter++)
285     mMultiplicities[tIter] = mRandom->Poisson(mAverageMultiplicities[tIter]);
286 }
287
288 void 
289 Event::DecayParticles()
290 {
291   Particle *tPart1=NULL, *tPart2=NULL, *tPart3=NULL, *tFather=NULL;
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     {
301       if (tFather->GetParticleType()->GetGamma() >= 0.0) {
302         if ((tFather->GetParticleType()->GetTable()) && (((DecayTable *) tFather->GetParticleType()->GetTable())->GetChannelCount()+1 > 0))
303           {
304             tDecayer->DecayParticle(tFather, &tPart1, &tPart2, &tPart3);
305 #ifndef _RESCALE_CHANNELS_
306             if ((!tPart1) || (!tPart2)) {
307               tCount++;
308               if (tPart1) delete tPart1;
309               if (tPart2) delete tPart2;
310               continue;
311             }
312 #endif
313 #ifdef _NO_THREE_BODY_DECAYS_
314             if (tPart3){
315               tCount++;
316               if (tPart1) delete tPart1;
317               if (tPart2) delete tPart2;
318               continue;
319             }
320 #endif
321 #ifdef _OMIT_TWO_BODY_
322             if (tPart1 && tPart2 && !tPart3){
323               tCount++;
324               delete tPart1;
325               delete tPart2;
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           }
349       }
350       tCount++;
351     }
352
353   delete tDecayer;
354 }
355
356 void 
357 Event::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
368 void           
369 Event::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   }
393   try {
394     mFOHSlocation = sRPInstance->getPar("FOHSLocation");
395   }
396   catch (STR e) {
397   }
398         
399 }
400