Fix Coverity
[u/mrichter/AliRoot.git] / TTherminator / TTherminator.cxx
CommitLineData
2e967919 1///////////////////////////////////////////////////////////////////////////////
2// //
3// TTherminator: global interface class to the Therminator model. //
4// Initialized global variables and runs the event generation //
5// //
6///////////////////////////////////////////////////////////////////////////////
7/******************************************************************************
8 * T H E R M I N A T O R *
9 * THERMal heavy-IoN generATOR *
10 * version 1.0 *
11 * *
12 * Authors of the model: Wojciech Broniowski, Wojciech.Broniowski@ifj.edu.pl, *
13 * Wojciech Florkowski, Wojciech.Florkowski@ifj.edu.pl *
14 * Authors of the code: Adam Kisiel, kisiel@if.pw.edu.pl *
15 * Tomasz Taluc, ttaluc@if.pw.edu.pl *
16 * Code designers: Adam Kisiel, Tomasz Taluc, Wojciech Broniowski, *
17 * Wojciech Florkowski *
18 * *
19 * For the detailed description of the program and furhter references *
20 * to the description of the model plesase refer to: nucl-th/0504047, *
21 * accessibile at: http://www.arxiv.org/nucl-th/0504047 *
22 * *
23 * Homepage: http://hirg.if.pw.edu.pl/en/therminator/ *
24 * *
25 * This code can be freely used and redistributed. However if you decide to *
26 * make modifications to the code, please contact the authors, especially *
27 * if you plan to publish the results obtained with such modified code. *
28 * Any publication of results obtained using this code must include the *
29 * reference to nucl-th/0504047 and the published version of it, when *
30 * available. *
31 * *
32 *****************************************************************************/
33#include <TTherminator.h>
34#include <fstream>
35#include <TDatabasePDG.h>
36#include <TParticle.h>
37#include <TClonesArray.h>
0ff6d0b8 38#include <TClass.h>
2e967919 39
40ReadPar *sRPInstance;
41STR sRPFileName;
42int sNumEvents;
43int sRunType;
44int sTables;
45int sModel;
46int sIntegrateSample;
47
48ClassImp(TTherminator)
49
50TTherminator::TTherminator():
51 fCalka(0),
52 fEvent(0),
53 fPartDB(0)
54{
55 // Default constructor
56 fPartDB = new ParticleDB();
57}
58TTherminator::TTherminator(const TTherminator & therm) :
c61a7285 59 TGenerator(therm),
ef6c3660 60 fCalka(new Integrator(*therm.fCalka)),
61 fEvent(new Event(*therm.fEvent)),
5c8cdc24 62 fPartDB(new ParticleDB())
2e967919 63{
64 // Copy constructor
5c8cdc24 65 // fPartDB = new ParticleDB();
66}
67TTherminator& TTherminator::operator=(const TTherminator & therm)
68{
69 if (this != &therm) {
70 fCalka = therm.fCalka;
71 fEvent = therm.fEvent;
72 delete fPartDB;
73 fPartDB = new ParticleDB();
74 }
75
76 return *this;
2e967919 77}
5c8cdc24 78
2e967919 79TTherminator::~TTherminator()
80{
81 // Destructor
82 if (sRPInstance) delete sRPInstance;
83 if (fEvent) delete fEvent;
84 if (fCalka) delete fCalka;
85 if (fPartDB) delete fPartDB;
86}
87
88void TTherminator::ReadParameters()
89{
90 // Read in global model parameters
91 STR tModel;
92 STR tTable;
93
94 try {
95 sNumEvents = atoi(sRPInstance->getPar("NumberOfEvents").Data());
96 tModel = sRPInstance->getPar("FreezeOutModel");
97 if (tModel.Contains("SingleFreezeOut"))
98 sModel = 0;
99/*MCH begin*/
100 else if (tModel.Contains("Lhyquid3D"))
101 sModel = 10;
102 else if (tModel.Contains("Lhyquid2D"))
103 sModel = 11;
104/*MCH end*/
105 else if (tModel.Contains("BlastWaveVTDelay"))
106 sModel = 6;
107 else if (tModel.Contains("BlastWaveVT"))
108 sModel = 2;
109 else if (tModel.Contains("BlastWaveVLinearFormation"))
110 sModel = 7;
111 else if (tModel.Contains("BlastWaveVLinearDelay"))
112 sModel = 8;
113 else if (tModel.Contains("BlastWaveVLinear"))
114 sModel = 4;
115 else if (tModel.Contains("FiniteHydro"))
116 sModel = 5;
117 else {
118 PRINT_MESSAGE("Unknown model type: " << tModel.Data());
119 PRINT_MESSAGE("Please provide the proper name");
120 exit(0);
121 }
122 if (atoi(sRPInstance->getPar("Randomize").Data()) == 1)
123 sRunType = 4;
124 else
125 sRunType = 3;
126 tTable = sRPInstance->getPar("TableType");
127 if (tTable.Contains("Mathematica"))
128 sTables = 0;
129 else if (tTable.Contains("SHARE"))
130 sTables = 1;
131 else {
132 PRINT_MESSAGE("Unknown table type: " << tTable.Data());
133 exit(0);
134 }
135 sIntegrateSample = atoi(sRPInstance->getPar("NumberOfIntegrateSamples").Data());
ae89e57a 136 fInputDir = sRPInstance->getPar("InputDirSHARE");
2e967919 137 }
138 catch (STR tError) {
139 PRINT_DEBUG_1("RunBFPW::ReadParameters - Caught exception " << tError);
140 PRINT_MESSAGE("Did not find one of the neccessary parameters in the parameters file.");
141 exit(0);
142 }
143}
144
145void TTherminator::Initialize(){
146 // Initialize global variables
147 // and particle and decay tables
148 Parser *tParser;
149
150 sRPFileName = "therminator.in";
151 sRPInstance = new ReadPar(sRPFileName.Data());
152
153 ReadParameters();
154
155 tParser = new Parser(fPartDB);
156
157 tParser->ReadShare();
158
159 fCalka = new Integrator(sIntegrateSample);
160 fCalka->ReadMultInteg(fPartDB);
161
162
163 // Read in particle table from share
164 // and add missing particle type to TDatabasePDG
165
166 char str[50];
167 char str1[200];
168
169 // ParticleType *tPartBuf;
170
171 TDatabasePDG *tInstance = TDatabasePDG::Instance();
172 TParticlePDG *tParticleType;
173
174 // AliWarning(Form("Reading particle types from particles.data"));
ae89e57a 175
176 ifstream in((fInputDir+"/"+"particles.data").Data());
177 // ifstream in("particles.data");
2e967919 178
179 int charge;
180
181 int number=0;
182 if ((in) && (in.is_open()))
183 {
184 //START OF HEAD-LINE
185 in.ignore(200,'\n');
186 in.ignore(200,'\n');
187 in.ignore(200,'\n');
188 //END OF HEAD-LINE
189
190 while (in>>str)
191 {
192 if (/*(*str == '#')||*/(*str<65)||(*str>122))
193 {
194 in.getline(str1,200);
195 continue;
196 }
197 double mass, gamma, spin, tI3, tI, q, s, aq, as, c, ac, mc;
198
199 in>>mass>>gamma>>spin>>tI>>tI3>>q>>s>>aq>>as>>c>>ac>>mc;
200 number++;
201 tParticleType = tInstance->GetParticle((int) mc);
202 if (!tParticleType) {
203 charge = 0;
204 if (strstr(str, "plu")) charge = 1;
205 if (strstr(str, "min")) charge = -1;
206 if (strstr(str, "plb")) charge = -1;
207 if (strstr(str, "mnb")) charge = 1;
208 if (strstr(str, "plp")) charge = 2;
209 if (strstr(str, "ppb")) charge = -2;
210 tInstance->AddParticle(str, str, mass, gamma == 0.0 ? 1:0, gamma, charge , "meson", (int) mc);
211 // printf("Added particle %s with PDG PID %d charge %d", str, (int) mc, charge);
212 //AliWarning(Form("Added particle %s with PDG PID %d charge %d", str, (int) mc, charge));
213 // AliWarning(Form("Quantum numbers q s c aq as ac tI3 %lf %lf %lf %lf %lf %lf %lf", q, s, c, aq, as, ac, tI3));
214
215 }
216 }
217 in.close();
218 }
09b7c66c 219
220 delete tParser;
2e967919 221}
222
223void TTherminator::GenerateEvent()
224{
225 // Run single event generation
226 if ((sRunType == 3) || (sRunType == 4))
227 {
228 if (sRunType==4)
229 fCalka->Randomize();
230
231 if (!fEvent) fEvent = new Event(fPartDB, fCalka);
232 if (sRunType == 4) fEvent->Randomize();
233
234 fEvent->GenerateEvent();
235 fEvent->DecayParticles();
236 // fEvent->WriteEvent(0);
237 }
238}
239
c61a7285 240Int_t TTherminator::ImportParticles(TClonesArray *particles, Option_t */*option*/)
2e967919 241{
242 // Import particles from a generated event into an external array
c9c8970c 243 const double kFmToGev = 0.197327;
244
2e967919 245 if (particles == 0) return 0;
246 TClonesArray &particlesR = *particles;
247 particlesR.Clear();
248 Int_t nump = 0;
249 if (!fEvent) return 0;
250 Int_t numpart = fEvent->GetParticleCount();
0ff6d0b8 251 printf("\n TTherminator: Therminator stack contains %d particles.\n", numpart);
2e967919 252 for (Int_t iPart=0; iPart<numpart; iPart++) {
253 Particle *tPart = fEvent->GetParticleOfCount(iPart);
254 Int_t tFather;
255 tFather = tPart->GetFather();
256
257 if (tFather> -1) {
258 TParticle *mother = (TParticle*) (particlesR.UncheckedAt(tFather+1));
259 mother->SetLastDaughter(iPart);
260 if (mother->GetFirstDaughter()==-1)
261 mother->SetFirstDaughter(iPart);
262 }
263 nump++;
264 // printf("Putting %d %d %lf %d\n", tPart->GetParticleType()->GetPDGCode(), iPart, tPart->px, tPart);
265 new (particlesR[iPart]) TParticle(tPart->GetParticleType()->GetPDGCode(), tPart->HadDecayed(),
0ff6d0b8 266 tFather, -1, -1, -1,
267 tPart->px, tPart->py, tPart->pz, tPart->GetEnergy() ,
268 tPart->rx*1.e-13*kFmToGev, tPart->ry*1.e-13*kFmToGev, tPart->rz*1.e-13*kFmToGev, tPart->rt*1.e-13*kFmToGev/3e10);
2e967919 269 particlesR[iPart]->SetUniqueID(iPart);
270 }
271
272 return nump;
273}
274
c61a7285 275TObjArray* TTherminator::ImportParticles(Option_t */*option*/)
2e967919 276{
277 // Import generated particles into an internal array
c9c8970c 278 const double kFmToGev = 0.197327;
279
2e967919 280 Int_t nump = 0;
281 fParticles->Clear();
282 if (!fEvent) return 0;
283 Int_t numpart = fEvent->GetParticleCount();
284 printf("\n TTherminator: Therminator stack contains %d particles.", numpart);
285 fEvent->InitParticleScan();
286 for (Int_t iPart=0; iPart<numpart; iPart++) {
287 Particle *tPart = fEvent->GetNextParticle();
288 Int_t tFather;
289 tFather = tPart->GetFather();
290
291 if (tFather> -1) {
292 TParticle *mother = (TParticle*) (fParticles->UncheckedAt(tFather));
293 mother->SetLastDaughter(iPart);
294 if (mother->GetFirstDaughter()==-1)
295 mother->SetFirstDaughter(iPart);
296 }
297 TParticle* p = new TParticle(tPart->GetParticleType()->GetPDGCode(), tPart->HadDecayed(),
298 tFather, -1, -1, -1,
299 tPart->px, tPart->py, tPart->pz, tPart->GetEnergy() ,
c9c8970c 300 tPart->rx*1.e-13*kFmToGev, tPart->ry*1.e-13*kFmToGev, tPart->rz*1.e-13*kFmToGev, tPart->rt*1.e-13*kFmToGev/3e10);
2e967919 301 p->SetUniqueID(iPart);
302 fParticles->Add(p);
303 }
304 nump = fEvent->GetParticleCount();
305
306 return fParticles;
307
308}