1 ///////////////////////////////////////////////////////////////////////////////
3 // TTherminator: global interface class to the Therminator model. //
4 // Initialized global variables and runs the event generation //
6 ///////////////////////////////////////////////////////////////////////////////
7 /******************************************************************************
8 * T H E R M I N A T O R *
9 * THERMal heavy-IoN generATOR *
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 *
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 *
23 * Homepage: http://hirg.if.pw.edu.pl/en/therminator/ *
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 *
32 *****************************************************************************/
33 #include <TTherminator.h>
35 #include <TDatabasePDG.h>
36 #include <TParticle.h>
37 #include <TClonesArray.h>
48 ClassImp(TTherminator)
50 TTherminator::TTherminator():
55 // Default constructor
56 fPartDB = new ParticleDB();
58 TTherminator::TTherminator(const TTherminator & therm) :
65 // fPartDB = new ParticleDB();
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();
73 TTherminator& TTherminator::operator=(const TTherminator & therm)
76 fCalka = therm.fCalka;
77 fEvent = therm.fEvent;
79 fPartDB = new ParticleDB();
85 TTherminator::~TTherminator()
88 if (sRPInstance) delete sRPInstance;
89 if (fEvent) delete fEvent;
90 if (fCalka) delete fCalka;
91 if (fPartDB) delete fPartDB;
94 void TTherminator::ReadParameters()
96 // Read in global model parameters
101 sNumEvents = atoi(sRPInstance->getPar("NumberOfEvents").Data());
102 tModel = sRPInstance->getPar("FreezeOutModel");
103 if (tModel.Contains("SingleFreezeOut"))
106 else if (tModel.Contains("Lhyquid3D"))
108 else if (tModel.Contains("Lhyquid2D"))
111 else if (tModel.Contains("BlastWaveVTDelay"))
113 else if (tModel.Contains("BlastWaveVT"))
115 else if (tModel.Contains("BlastWaveVLinearFormation"))
117 else if (tModel.Contains("BlastWaveVLinearDelay"))
119 else if (tModel.Contains("BlastWaveVLinear"))
121 else if (tModel.Contains("FiniteHydro"))
124 PRINT_MESSAGE("Unknown model type: " << tModel.Data());
125 PRINT_MESSAGE("Please provide the proper name");
128 if (atoi(sRPInstance->getPar("Randomize").Data()) == 1)
132 tTable = sRPInstance->getPar("TableType");
133 if (tTable.Contains("Mathematica"))
135 else if (tTable.Contains("SHARE"))
138 PRINT_MESSAGE("Unknown table type: " << tTable.Data());
141 sIntegrateSample = atoi(sRPInstance->getPar("NumberOfIntegrateSamples").Data());
142 fInputDir = sRPInstance->getPar("InputDirSHARE");
145 PRINT_DEBUG_1("RunBFPW::ReadParameters - Caught exception " << tError);
146 PRINT_MESSAGE("Did not find one of the neccessary parameters in the parameters file.");
151 void TTherminator::Initialize(){
152 // Initialize global variables
153 // and particle and decay tables
156 sRPFileName = "therminator.in";
157 sRPInstance = new ReadPar(sRPFileName.Data());
161 tParser = new Parser(fPartDB);
163 tParser->ReadShare();
165 fCalka = new Integrator(sIntegrateSample);
166 fCalka->ReadMultInteg(fPartDB);
169 // Read in particle table from share
170 // and add missing particle type to TDatabasePDG
175 // ParticleType *tPartBuf;
177 TDatabasePDG *tInstance = TDatabasePDG::Instance();
178 TParticlePDG *tParticleType;
180 // AliWarning(Form("Reading particle types from particles.data"));
182 ifstream in((fInputDir+"/"+"particles.data").Data());
183 // ifstream in("particles.data");
188 if ((in) && (in.is_open()))
198 if (/*(*str == '#')||*/(*str<65)||(*str>122))
200 in.getline(str1,200);
203 double mass, gamma, spin, tI3, tI, q, s, aq, as, c, ac, mc;
205 in>>mass>>gamma>>spin>>tI>>tI3>>q>>s>>aq>>as>>c>>ac>>mc;
207 tParticleType = tInstance->GetParticle((int) mc);
208 if (!tParticleType) {
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));
229 void TTherminator::GenerateEvent()
231 // Run single event generation
232 if ((sRunType == 3) || (sRunType == 4))
237 if (!fEvent) fEvent = new Event(fPartDB, fCalka);
238 if (sRunType == 4) fEvent->Randomize();
240 fEvent->GenerateEvent();
241 fEvent->DecayParticles();
242 // fEvent->WriteEvent(0);
246 Int_t TTherminator::ImportParticles(TClonesArray *particles, Option_t */*option*/)
248 // Import particles from a generated event into an external array
249 const double kFmToGev = 0.197327;
251 if (particles == 0) return 0;
252 TClonesArray &particlesR = *particles;
255 if (!fEvent) return 0;
256 Int_t numpart = fEvent->GetParticleCount();
257 printf("\n TTherminator: Therminator stack contains %d particles.\n", numpart);
258 for (Int_t iPart=0; iPart<numpart; iPart++) {
259 Particle *tPart = fEvent->GetParticleOfCount(iPart);
261 tFather = tPart->GetFather();
264 TParticle *mother = (TParticle*) (particlesR.UncheckedAt(tFather+1));
265 mother->SetLastDaughter(iPart);
266 if (mother->GetFirstDaughter()==-1)
267 mother->SetFirstDaughter(iPart);
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(),
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);
275 particlesR[iPart]->SetUniqueID(iPart);
281 TObjArray* TTherminator::ImportParticles(Option_t */*option*/)
283 // Import generated particles into an internal array
284 const double kFmToGev = 0.197327;
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();
295 tFather = tPart->GetFather();
298 TParticle *mother = (TParticle*) (fParticles->UncheckedAt(tFather));
299 mother->SetLastDaughter(iPart);
300 if (mother->GetFirstDaughter()==-1)
301 mother->SetFirstDaughter(iPart);
303 TParticle* p = new TParticle(tPart->GetParticleType()->GetPDGCode(), tPart->HadDecayed(),
305 tPart->px, tPart->py, tPart->pz, tPart->GetEnergy() ,
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);
307 p->SetUniqueID(iPart);
310 nump = fEvent->GetParticleCount();