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();
67 TTherminator::~TTherminator()
70 if (sRPInstance) delete sRPInstance;
71 if (fEvent) delete fEvent;
72 if (fCalka) delete fCalka;
73 if (fPartDB) delete fPartDB;
76 void TTherminator::ReadParameters()
78 // Read in global model parameters
83 sNumEvents = atoi(sRPInstance->getPar("NumberOfEvents").Data());
84 tModel = sRPInstance->getPar("FreezeOutModel");
85 if (tModel.Contains("SingleFreezeOut"))
88 else if (tModel.Contains("Lhyquid3D"))
90 else if (tModel.Contains("Lhyquid2D"))
93 else if (tModel.Contains("BlastWaveVTDelay"))
95 else if (tModel.Contains("BlastWaveVT"))
97 else if (tModel.Contains("BlastWaveVLinearFormation"))
99 else if (tModel.Contains("BlastWaveVLinearDelay"))
101 else if (tModel.Contains("BlastWaveVLinear"))
103 else if (tModel.Contains("FiniteHydro"))
106 PRINT_MESSAGE("Unknown model type: " << tModel.Data());
107 PRINT_MESSAGE("Please provide the proper name");
110 if (atoi(sRPInstance->getPar("Randomize").Data()) == 1)
114 tTable = sRPInstance->getPar("TableType");
115 if (tTable.Contains("Mathematica"))
117 else if (tTable.Contains("SHARE"))
120 PRINT_MESSAGE("Unknown table type: " << tTable.Data());
123 sIntegrateSample = atoi(sRPInstance->getPar("NumberOfIntegrateSamples").Data());
124 fInputDir = sRPInstance->getPar("InputDirSHARE");
127 PRINT_DEBUG_1("RunBFPW::ReadParameters - Caught exception " << tError);
128 PRINT_MESSAGE("Did not find one of the neccessary parameters in the parameters file.");
133 void TTherminator::Initialize(){
134 // Initialize global variables
135 // and particle and decay tables
138 sRPFileName = "therminator.in";
139 sRPInstance = new ReadPar(sRPFileName.Data());
143 tParser = new Parser(fPartDB);
145 tParser->ReadShare();
147 fCalka = new Integrator(sIntegrateSample);
148 fCalka->ReadMultInteg(fPartDB);
151 // Read in particle table from share
152 // and add missing particle type to TDatabasePDG
157 // ParticleType *tPartBuf;
159 TDatabasePDG *tInstance = TDatabasePDG::Instance();
160 TParticlePDG *tParticleType;
162 // AliWarning(Form("Reading particle types from particles.data"));
164 ifstream in((fInputDir+"/"+"particles.data").Data());
165 // ifstream in("particles.data");
170 if ((in) && (in.is_open()))
180 if (/*(*str == '#')||*/(*str<65)||(*str>122))
182 in.getline(str1,200);
185 double mass, gamma, spin, tI3, tI, q, s, aq, as, c, ac, mc;
187 in>>mass>>gamma>>spin>>tI>>tI3>>q>>s>>aq>>as>>c>>ac>>mc;
189 tParticleType = tInstance->GetParticle((int) mc);
190 if (!tParticleType) {
192 if (strstr(str, "plu")) charge = 1;
193 if (strstr(str, "min")) charge = -1;
194 if (strstr(str, "plb")) charge = -1;
195 if (strstr(str, "mnb")) charge = 1;
196 if (strstr(str, "plp")) charge = 2;
197 if (strstr(str, "ppb")) charge = -2;
198 tInstance->AddParticle(str, str, mass, gamma == 0.0 ? 1:0, gamma, charge , "meson", (int) mc);
199 // printf("Added particle %s with PDG PID %d charge %d", str, (int) mc, charge);
200 //AliWarning(Form("Added particle %s with PDG PID %d charge %d", str, (int) mc, charge));
201 // 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));
209 void TTherminator::GenerateEvent()
211 // Run single event generation
212 if ((sRunType == 3) || (sRunType == 4))
217 if (!fEvent) fEvent = new Event(fPartDB, fCalka);
218 if (sRunType == 4) fEvent->Randomize();
220 fEvent->GenerateEvent();
221 fEvent->DecayParticles();
222 // fEvent->WriteEvent(0);
226 Int_t TTherminator::ImportParticles(TClonesArray *particles, Option_t */*option*/)
228 // Import particles from a generated event into an external array
229 const double kFmToGev = 0.197327;
231 if (particles == 0) return 0;
232 TClonesArray &particlesR = *particles;
235 if (!fEvent) return 0;
236 Int_t numpart = fEvent->GetParticleCount();
237 printf("\n TTherminator: Therminator stack contains %d particles.\n", numpart);
238 for (Int_t iPart=0; iPart<numpart; iPart++) {
239 Particle *tPart = fEvent->GetParticleOfCount(iPart);
241 tFather = tPart->GetFather();
244 TParticle *mother = (TParticle*) (particlesR.UncheckedAt(tFather+1));
245 mother->SetLastDaughter(iPart);
246 if (mother->GetFirstDaughter()==-1)
247 mother->SetFirstDaughter(iPart);
250 // printf("Putting %d %d %lf %d\n", tPart->GetParticleType()->GetPDGCode(), iPart, tPart->px, tPart);
251 new (particlesR[iPart]) TParticle(tPart->GetParticleType()->GetPDGCode(), tPart->HadDecayed(),
253 tPart->px, tPart->py, tPart->pz, tPart->GetEnergy() ,
254 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);
255 particlesR[iPart]->SetUniqueID(iPart);
261 TObjArray* TTherminator::ImportParticles(Option_t */*option*/)
263 // Import generated particles into an internal array
264 const double kFmToGev = 0.197327;
268 if (!fEvent) return 0;
269 Int_t numpart = fEvent->GetParticleCount();
270 printf("\n TTherminator: Therminator stack contains %d particles.", numpart);
271 fEvent->InitParticleScan();
272 for (Int_t iPart=0; iPart<numpart; iPart++) {
273 Particle *tPart = fEvent->GetNextParticle();
275 tFather = tPart->GetFather();
278 TParticle *mother = (TParticle*) (fParticles->UncheckedAt(tFather));
279 mother->SetLastDaughter(iPart);
280 if (mother->GetFirstDaughter()==-1)
281 mother->SetFirstDaughter(iPart);
283 TParticle* p = new TParticle(tPart->GetParticleType()->GetPDGCode(), tPart->HadDecayed(),
285 tPart->px, tPart->py, tPart->pz, tPart->GetEnergy() ,
286 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);
287 p->SetUniqueID(iPart);
290 nump = fEvent->GetParticleCount();