]>
Commit | Line | Data |
---|---|---|
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 | ||
39 | ReadPar *sRPInstance; | |
40 | STR sRPFileName; | |
41 | int sNumEvents; | |
42 | int sRunType; | |
43 | int sTables; | |
44 | int sModel; | |
45 | int sIntegrateSample; | |
46 | ||
47 | ClassImp(TTherminator) | |
48 | ||
49 | TTherminator::TTherminator(): | |
50 | fCalka(0), | |
51 | fEvent(0), | |
52 | fPartDB(0) | |
53 | { | |
54 | // Default constructor | |
55 | fPartDB = new ParticleDB(); | |
56 | } | |
57 | TTherminator::TTherminator(const TTherminator & therm) : | |
58 | TGenerator(therm), | |
59 | fCalka(0), | |
60 | fEvent(0), | |
61 | fPartDB(0) | |
62 | { | |
63 | // Copy constructor | |
64 | fPartDB = new ParticleDB(); | |
65 | } | |
66 | TTherminator::~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 | ||
75 | void 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()); | |
123 | fInputDir = sRPInstance->getPar("InputDirSHARE"); | |
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 | ||
132 | void 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")); | |
162 | ||
163 | ifstream in((fInputDir+"/"+"particles.data").Data()); | |
164 | // ifstream in("particles.data"); | |
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 | ||
208 | void 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 | ||
225 | Int_t TTherminator::ImportParticles(TClonesArray *particles, Option_t */*option*/) | |
226 | { | |
227 | // Import particles from a generated event into an external array | |
228 | const double kFmToGev = 0.197327; | |
229 | ||
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() , | |
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); | |
254 | particlesR[iPart]->SetUniqueID(iPart); | |
255 | } | |
256 | ||
257 | return nump; | |
258 | } | |
259 | ||
260 | TObjArray* TTherminator::ImportParticles(Option_t */*option*/) | |
261 | { | |
262 | // Import generated particles into an internal array | |
263 | const double kFmToGev = 0.197327; | |
264 | ||
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() , | |
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); | |
286 | p->SetUniqueID(iPart); | |
287 | fParticles->Add(p); | |
288 | } | |
289 | nump = fEvent->GetParticleCount(); | |
290 | ||
291 | return fParticles; | |
292 | ||
293 | } |