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