Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[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) :
1aed17af 59 TGenerator(therm),
60 fCalka(0),
61 fEvent(0),
62 fPartDB(0)
2e967919 63{
64 // Copy constructor
5c8cdc24 65 // fPartDB = new ParticleDB();
1aed17af 66 if (fCalka) delete fCalka;
67 fCalka = new Integrator(*therm.fCalka);
68 if (fEvent) delete fEvent;
69 fEvent = new Event(*therm.fEvent);
70 if (fPartDB) delete fPartDB;
71 fPartDB = new ParticleDB();
5c8cdc24 72}
73TTherminator& TTherminator::operator=(const TTherminator & therm)
74{
75 if (this != &therm) {
76 fCalka = therm.fCalka;
77 fEvent = therm.fEvent;
78 delete fPartDB;
79 fPartDB = new ParticleDB();
80 }
81
82 return *this;
2e967919 83}
5c8cdc24 84
2e967919 85TTherminator::~TTherminator()
86{
87 // Destructor
88 if (sRPInstance) delete sRPInstance;
89 if (fEvent) delete fEvent;
90 if (fCalka) delete fCalka;
91 if (fPartDB) delete fPartDB;
92}
93
94void TTherminator::ReadParameters()
95{
96 // Read in global model parameters
97 STR tModel;
98 STR tTable;
99
100 try {
101 sNumEvents = atoi(sRPInstance->getPar("NumberOfEvents").Data());
102 tModel = sRPInstance->getPar("FreezeOutModel");
103 if (tModel.Contains("SingleFreezeOut"))
104 sModel = 0;
105/*MCH begin*/
106 else if (tModel.Contains("Lhyquid3D"))
107 sModel = 10;
108 else if (tModel.Contains("Lhyquid2D"))
109 sModel = 11;
110/*MCH end*/
111 else if (tModel.Contains("BlastWaveVTDelay"))
112 sModel = 6;
113 else if (tModel.Contains("BlastWaveVT"))
114 sModel = 2;
115 else if (tModel.Contains("BlastWaveVLinearFormation"))
116 sModel = 7;
117 else if (tModel.Contains("BlastWaveVLinearDelay"))
118 sModel = 8;
119 else if (tModel.Contains("BlastWaveVLinear"))
120 sModel = 4;
121 else if (tModel.Contains("FiniteHydro"))
122 sModel = 5;
123 else {
124 PRINT_MESSAGE("Unknown model type: " << tModel.Data());
125 PRINT_MESSAGE("Please provide the proper name");
126 exit(0);
127 }
128 if (atoi(sRPInstance->getPar("Randomize").Data()) == 1)
129 sRunType = 4;
130 else
131 sRunType = 3;
132 tTable = sRPInstance->getPar("TableType");
133 if (tTable.Contains("Mathematica"))
134 sTables = 0;
135 else if (tTable.Contains("SHARE"))
136 sTables = 1;
137 else {
138 PRINT_MESSAGE("Unknown table type: " << tTable.Data());
139 exit(0);
140 }
141 sIntegrateSample = atoi(sRPInstance->getPar("NumberOfIntegrateSamples").Data());
ae89e57a 142 fInputDir = sRPInstance->getPar("InputDirSHARE");
2e967919 143 }
144 catch (STR tError) {
145 PRINT_DEBUG_1("RunBFPW::ReadParameters - Caught exception " << tError);
146 PRINT_MESSAGE("Did not find one of the neccessary parameters in the parameters file.");
147 exit(0);
148 }
149}
150
151void TTherminator::Initialize(){
152 // Initialize global variables
153 // and particle and decay tables
154 Parser *tParser;
155
156 sRPFileName = "therminator.in";
157 sRPInstance = new ReadPar(sRPFileName.Data());
158
159 ReadParameters();
160
161 tParser = new Parser(fPartDB);
162
163 tParser->ReadShare();
164
165 fCalka = new Integrator(sIntegrateSample);
166 fCalka->ReadMultInteg(fPartDB);
167
168
169 // Read in particle table from share
170 // and add missing particle type to TDatabasePDG
171
172 char str[50];
173 char str1[200];
174
175 // ParticleType *tPartBuf;
176
177 TDatabasePDG *tInstance = TDatabasePDG::Instance();
178 TParticlePDG *tParticleType;
179
180 // AliWarning(Form("Reading particle types from particles.data"));
ae89e57a 181
182 ifstream in((fInputDir+"/"+"particles.data").Data());
183 // ifstream in("particles.data");
2e967919 184
185 int charge;
186
187 int number=0;
188 if ((in) && (in.is_open()))
189 {
190 //START OF HEAD-LINE
191 in.ignore(200,'\n');
192 in.ignore(200,'\n');
193 in.ignore(200,'\n');
194 //END OF HEAD-LINE
195
196 while (in>>str)
197 {
198 if (/*(*str == '#')||*/(*str<65)||(*str>122))
199 {
200 in.getline(str1,200);
201 continue;
202 }
203 double mass, gamma, spin, tI3, tI, q, s, aq, as, c, ac, mc;
204
205 in>>mass>>gamma>>spin>>tI>>tI3>>q>>s>>aq>>as>>c>>ac>>mc;
206 number++;
207 tParticleType = tInstance->GetParticle((int) mc);
208 if (!tParticleType) {
209 charge = 0;
210 if (strstr(str, "plu")) charge = 1;
211 if (strstr(str, "min")) charge = -1;
212 if (strstr(str, "plb")) charge = -1;
213 if (strstr(str, "mnb")) charge = 1;
214 if (strstr(str, "plp")) charge = 2;
215 if (strstr(str, "ppb")) charge = -2;
216 tInstance->AddParticle(str, str, mass, gamma == 0.0 ? 1:0, gamma, charge , "meson", (int) mc);
217 // printf("Added particle %s with PDG PID %d charge %d", str, (int) mc, charge);
218 //AliWarning(Form("Added particle %s with PDG PID %d charge %d", str, (int) mc, charge));
219 // 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));
220
221 }
222 }
223 in.close();
224 }
09b7c66c 225
226 delete tParser;
2e967919 227}
228
229void TTherminator::GenerateEvent()
230{
231 // Run single event generation
232 if ((sRunType == 3) || (sRunType == 4))
233 {
234 if (sRunType==4)
235 fCalka->Randomize();
236
237 if (!fEvent) fEvent = new Event(fPartDB, fCalka);
238 if (sRunType == 4) fEvent->Randomize();
239
240 fEvent->GenerateEvent();
241 fEvent->DecayParticles();
242 // fEvent->WriteEvent(0);
243 }
244}
245
c61a7285 246Int_t TTherminator::ImportParticles(TClonesArray *particles, Option_t */*option*/)
2e967919 247{
248 // Import particles from a generated event into an external array
c9c8970c 249 const double kFmToGev = 0.197327;
250
2e967919 251 if (particles == 0) return 0;
252 TClonesArray &particlesR = *particles;
253 particlesR.Clear();
254 Int_t nump = 0;
255 if (!fEvent) return 0;
256 Int_t numpart = fEvent->GetParticleCount();
0ff6d0b8 257 printf("\n TTherminator: Therminator stack contains %d particles.\n", numpart);
2e967919 258 for (Int_t iPart=0; iPart<numpart; iPart++) {
259 Particle *tPart = fEvent->GetParticleOfCount(iPart);
260 Int_t tFather;
261 tFather = tPart->GetFather();
262
263 if (tFather> -1) {
264 TParticle *mother = (TParticle*) (particlesR.UncheckedAt(tFather+1));
265 mother->SetLastDaughter(iPart);
266 if (mother->GetFirstDaughter()==-1)
267 mother->SetFirstDaughter(iPart);
268 }
269 nump++;
270 // printf("Putting %d %d %lf %d\n", tPart->GetParticleType()->GetPDGCode(), iPart, tPart->px, tPart);
271 new (particlesR[iPart]) TParticle(tPart->GetParticleType()->GetPDGCode(), tPart->HadDecayed(),
0ff6d0b8 272 tFather, -1, -1, -1,
273 tPart->px, tPart->py, tPart->pz, tPart->GetEnergy() ,
274 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 275 particlesR[iPart]->SetUniqueID(iPart);
276 }
277
278 return nump;
279}
280
c61a7285 281TObjArray* TTherminator::ImportParticles(Option_t */*option*/)
2e967919 282{
283 // Import generated particles into an internal array
c9c8970c 284 const double kFmToGev = 0.197327;
285
2e967919 286 Int_t nump = 0;
287 fParticles->Clear();
288 if (!fEvent) return 0;
289 Int_t numpart = fEvent->GetParticleCount();
290 printf("\n TTherminator: Therminator stack contains %d particles.", numpart);
291 fEvent->InitParticleScan();
292 for (Int_t iPart=0; iPart<numpart; iPart++) {
293 Particle *tPart = fEvent->GetNextParticle();
294 Int_t tFather;
295 tFather = tPart->GetFather();
296
297 if (tFather> -1) {
298 TParticle *mother = (TParticle*) (fParticles->UncheckedAt(tFather));
299 mother->SetLastDaughter(iPart);
300 if (mother->GetFirstDaughter()==-1)
301 mother->SetFirstDaughter(iPart);
302 }
303 TParticle* p = new TParticle(tPart->GetParticleType()->GetPDGCode(), tPart->HadDecayed(),
304 tFather, -1, -1, -1,
305 tPart->px, tPart->py, tPart->pz, tPart->GetEnergy() ,
c9c8970c 306 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 307 p->SetUniqueID(iPart);
308 fParticles->Add(p);
309 }
310 nump = fEvent->GetParticleCount();
311
312 return fParticles;
313
314}