]>
Commit | Line | Data |
---|---|---|
2e967919 | 1 | // ALICE event generator based on the THERMINATOR model |
2 | // It reads the test output of the model and puts it onto | |
3 | // the stack | |
ae89e57a | 4 | // It has an option to use the Lhyquid3D input freeze-out |
5 | // hypersurface | |
2e967919 | 6 | // Author: Adam.Kisiel@cern.ch |
7 | ||
8 | #include <iostream> | |
9 | #include <fstream> | |
10 | #include <sstream> | |
11 | #include <TClonesArray.h> | |
12 | #include <TMCProcess.h> | |
13 | #include <TDatabasePDG.h> | |
14 | #include <TParticle.h> | |
15 | ||
16 | #include "AliConst.h" | |
17 | #include "AliDecayer.h" | |
18 | #include "AliGenEventHeader.h" | |
9abb118d | 19 | #include "AliGenHijingEventHeader.h" |
2e967919 | 20 | #include "AliGenTherminator.h" |
21 | #include "AliLog.h" | |
22 | #include "AliRun.h" | |
23 | ||
24 | ClassImp(AliGenTherminator) | |
25 | ||
26 | using namespace std; | |
27 | ||
28 | AliGenTherminator::AliGenTherminator(): | |
29 | AliGenMC(), | |
30 | fNt(0), | |
31 | fEventNumber(0), | |
32 | fFileName(""), | |
33 | fFreezeOutModel(""), | |
34 | fTemperature(0.1656), | |
35 | fMiuI(-0.0009), | |
36 | fMiuS(0.0), | |
37 | fMiuB(0.0008), | |
38 | fAlfaRange(6.0), | |
39 | fRapRange(4.0), | |
40 | fRhoMax(8.0), | |
41 | fTau(8.0), | |
42 | fBWA(0.0), | |
43 | fBWVt(1.41), | |
44 | fBWDelay(0.0) | |
45 | { | |
46 | // Default constructor | |
2e967919 | 47 | } |
48 | AliGenTherminator::AliGenTherminator(Int_t npart): | |
49 | AliGenMC(npart), | |
50 | fNt(0), | |
51 | fEventNumber(0), | |
52 | fFileName(""), | |
53 | fFreezeOutModel(""), | |
54 | fTemperature(0.1656), | |
55 | fMiuI(-0.0009), | |
56 | fMiuS(0.0), | |
57 | fMiuB(0.0008), | |
58 | fAlfaRange(6.0), | |
59 | fRapRange(4.0), | |
60 | fRhoMax(8.0), | |
61 | fTau(8.0), | |
62 | fBWA(0.0), | |
63 | fBWVt(1.41), | |
64 | fBWDelay(0.0) | |
65 | { | |
66 | // Constructor specifying the size of the particle table | |
2e967919 | 67 | fNprimaries = 0; |
68 | } | |
69 | ||
70 | AliGenTherminator::~AliGenTherminator() | |
71 | { | |
72 | // AliGenMC::~AliGenMC(); | |
73 | // if (fTherminator) delete fTherminator; | |
74 | } | |
75 | ||
76 | void AliGenTherminator::Generate() | |
77 | { | |
78 | // Run single event generation with the Therminator model | |
79 | AliWarning("Generating event from AliGenTherminator"); | |
80 | ||
81 | Float_t polar[3] = {0,0,0}; | |
82 | Float_t origin[3] = {0,0,0}; | |
83 | Float_t origin0[3] = {0,0,0}; | |
84 | Float_t p[3]; | |
85 | Float_t mass, energy; | |
86 | ||
87 | Int_t nt = 0; | |
88 | Int_t j, kf, ks, imo; | |
89 | kf = 0; | |
90 | ||
91 | Vertex(); | |
92 | for (j=0; j < 3; j++) origin0[j] = fVertex[j]; | |
93 | ||
94 | // Generate one event | |
95 | ||
96 | ((TTherminator *) fMCEvGen)->GenerateEvent(); | |
97 | AliWarning("Generated"); | |
8507138f | 98 | ((TTherminator *) fMCEvGen)->ImportParticles(&fParticles); |
2e967919 | 99 | |
8507138f | 100 | Int_t np = fParticles.GetEntriesFast(); |
2e967919 | 101 | AliWarning(Form("Imported %d particles", np)); |
102 | ||
deb008f3 | 103 | Int_t *idsOnStack; |
104 | idsOnStack = new Int_t[np]; | |
105 | ||
2e967919 | 106 | TParticle *iparticle; |
9abb118d | 107 | Double_t evrot = gRandom->Rndm()*TMath::Pi(); |
2e967919 | 108 | |
109 | for (int i = 0; i < np; i++) { | |
8507138f | 110 | iparticle = (TParticle *) fParticles.At(i); |
2e967919 | 111 | Bool_t hasMother = (iparticle->GetFirstMother() >=0); |
112 | Bool_t hasDaughter = (iparticle->GetFirstDaughter() >=0); | |
113 | ||
c9a371fc | 114 | if (hasDaughter) { |
115 | // This particle has decayed | |
116 | // It will not be tracked | |
117 | // Add it only once with coorduinates not | |
118 | // smeared with primary vertex position | |
119 | ||
120 | kf = iparticle->GetPdgCode(); | |
121 | ks = iparticle->GetStatusCode(); | |
122 | Double_t aphi = TMath::ATan2(iparticle->Py(), iparticle->Px()); | |
123 | Double_t arho = TMath::Hypot(iparticle->Px(), iparticle->Py()); | |
124 | p[0] = arho*TMath::Cos(aphi + evrot); | |
125 | p[1] = arho*TMath::Sin(aphi + evrot); | |
126 | // p[0] = iparticle->Px(); | |
127 | // p[1] = iparticle->Py(); | |
128 | p[2] = iparticle->Pz(); | |
129 | mass = TDatabasePDG::Instance()->GetParticle(kf)->Mass(); | |
130 | energy = sqrt(mass*mass + p[0]*p[0] + p[1]*p[1] + p[2]*p[2]); | |
131 | ||
132 | Double_t vphi = TMath::ATan2(iparticle->Vy(), iparticle->Vx()); | |
133 | Double_t vrho = TMath::Hypot(iparticle->Vx(), iparticle->Vy()); | |
c9c8970c | 134 | origin[0] = vrho*TMath::Cos(vphi + evrot); |
135 | origin[1] = vrho*TMath::Sin(vphi + evrot); | |
136 | origin[2] = iparticle->Vz(); | |
c9a371fc | 137 | |
138 | imo = -1; | |
139 | TParticle* mother = 0; | |
140 | if (hasMother) { | |
141 | imo = iparticle->GetFirstMother(); | |
c9c8970c | 142 | mother = (TParticle *) fParticles.At(imo); |
c9a371fc | 143 | } // if has mother |
144 | Bool_t tFlag = (!hasDaughter); | |
145 | ||
d75811a6 | 146 | printf("Pushing Track %d with status %d mother %d\n", kf, tFlag, imo>=0?idsOnStack[imo]:imo); |
147 | PushTrack(tFlag,imo>=0?idsOnStack[imo]:imo,kf, | |
c9a371fc | 148 | p[0],p[1],p[2],energy, |
c9c8970c | 149 | origin[0],origin[1],origin[2],iparticle->T(), |
c9a371fc | 150 | polar[0],polar[1],polar[2], |
151 | hasMother ? kPDecay:kPNoProcess,nt); | |
deb008f3 | 152 | idsOnStack[i] = nt; |
c9a371fc | 153 | fNprimaries++; |
154 | KeepTrack(nt); | |
155 | } | |
156 | else { | |
157 | // This is a final state particle | |
158 | // It will be tracked | |
159 | // Add it TWICE to the stack !!! | |
160 | // First time with event-wide coordicates (for femtoscopy) - | |
161 | // this one will not be tracked | |
162 | // Second time with event-wide ccordiantes and vertex smearing | |
163 | // this one will be tracked | |
164 | ||
165 | kf = iparticle->GetPdgCode(); | |
166 | ks = iparticle->GetStatusCode(); | |
167 | Double_t aphi = TMath::ATan2(iparticle->Py(), iparticle->Px()); | |
168 | Double_t arho = TMath::Hypot(iparticle->Px(), iparticle->Py()); | |
169 | p[0] = arho*TMath::Cos(aphi + evrot); | |
170 | p[1] = arho*TMath::Sin(aphi + evrot); | |
171 | // p[0] = iparticle->Px(); | |
172 | // p[1] = iparticle->Py(); | |
173 | p[2] = iparticle->Pz(); | |
174 | mass = TDatabasePDG::Instance()->GetParticle(kf)->Mass(); | |
175 | energy = sqrt(mass*mass + p[0]*p[0] + p[1]*p[1] + p[2]*p[2]); | |
176 | ||
177 | Double_t vphi = TMath::ATan2(iparticle->Vy(), iparticle->Vx()); | |
178 | Double_t vrho = TMath::Hypot(iparticle->Vx(), iparticle->Vy()); | |
179 | origin[0] = vrho*TMath::Cos(vphi + evrot); | |
180 | origin[1] = vrho*TMath::Sin(vphi + evrot); | |
181 | origin[2] = iparticle->Vz(); | |
182 | ||
183 | imo = -1; | |
184 | TParticle* mother = 0; | |
185 | if (hasMother) { | |
186 | imo = iparticle->GetFirstMother(); | |
c9c8970c | 187 | mother = (TParticle *) fParticles.At(imo); |
c9a371fc | 188 | } // if has mother |
189 | Bool_t tFlag = (hasDaughter); | |
190 | ||
d75811a6 | 191 | printf("Found mother %i with true id %i\n", imo, imo>=0?idsOnStack[imo]:imo); |
192 | printf("Pushing Track %d with status %d mother %d\n", kf, tFlag, imo>=0?idsOnStack[imo]:imo); | |
193 | PushTrack(tFlag,imo>=0?idsOnStack[imo]:imo,kf, | |
c9a371fc | 194 | p[0],p[1],p[2],energy, |
c9c8970c | 195 | origin[0],origin[1],origin[2],iparticle->T(), |
c9a371fc | 196 | polar[0],polar[1],polar[2], |
197 | hasMother ? kPDecay:kPNoProcess,nt); | |
deb008f3 | 198 | idsOnStack[i] = nt; |
c9a371fc | 199 | fNprimaries++; |
200 | KeepTrack(nt); | |
201 | ||
202 | origin[0] = origin0[0]+vrho*TMath::Cos(vphi + evrot); | |
203 | origin[1] = origin0[1]+vrho*TMath::Sin(vphi + evrot); | |
204 | origin[2] = origin0[2]+iparticle->Vz(); | |
205 | ||
206 | imo = nt; | |
deb008f3 | 207 | // mother = (TParticle *) fParticles.At(nt); |
c9a371fc | 208 | tFlag = (!hasDaughter); |
209 | ||
210 | printf("Pushing Track %d with status %d mother %d\n", kf, tFlag, imo); | |
211 | PushTrack(tFlag,imo,kf, | |
212 | p[0],p[1],p[2],energy, | |
c9c8970c | 213 | origin[0],origin[1],origin[2],iparticle->T(), |
c9a371fc | 214 | polar[0],polar[1],polar[2], |
215 | hasMother ? kPDecay:kPNoProcess,nt); | |
216 | fNprimaries++; | |
217 | KeepTrack(nt); | |
218 | } | |
2e967919 | 219 | } |
220 | ||
deb008f3 | 221 | |
222 | ||
2e967919 | 223 | SetHighWaterMark(fNprimaries); |
224 | ||
225 | TArrayF eventVertex; | |
226 | eventVertex.Set(3); | |
227 | eventVertex[0] = origin0[0]; | |
228 | eventVertex[1] = origin0[1]; | |
229 | eventVertex[2] = origin0[2]; | |
230 | ||
9abb118d | 231 | // Builds the event header, to be called after each event |
232 | AliGenEventHeader* header = new AliGenHijingEventHeader("Therminator"); | |
233 | ||
2e967919 | 234 | // Header |
9abb118d | 235 | // AliGenEventHeader* header = new AliGenEventHeader("Therminator"); |
2e967919 | 236 | // Event Vertex |
9abb118d | 237 | // header->SetPrimaryVertex(eventVertex); |
238 | // header->SetNProduced(fNprimaries); | |
239 | ||
240 | ((AliGenHijingEventHeader*) header)->SetNProduced(fNprimaries); | |
241 | ((AliGenHijingEventHeader*) header)->SetPrimaryVertex(eventVertex); | |
242 | ((AliGenHijingEventHeader*) header)->SetImpactParameter(0.0); | |
243 | ((AliGenHijingEventHeader*) header)->SetTotalEnergy(0.0); | |
244 | ((AliGenHijingEventHeader*) header)->SetHardScatters(0); | |
245 | ((AliGenHijingEventHeader*) header)->SetParticipants(0, 0); | |
246 | ((AliGenHijingEventHeader*) header)->SetCollisions(0, 0, 0, 0); | |
247 | ((AliGenHijingEventHeader*) header)->SetSpectators(0, 0, 0, 0); | |
248 | ((AliGenHijingEventHeader*) header)->SetReactionPlaneAngle(evrot); | |
249 | ||
250 | ||
251 | // 4-momentum vectors of the triggered jets. | |
252 | // | |
253 | // Before final state gluon radiation. | |
254 | // TLorentzVector* jet1 = new TLorentzVector(fHijing->GetHINT1(21), | |
255 | // fHijing->GetHINT1(22), | |
256 | // fHijing->GetHINT1(23), | |
257 | // fHijing->GetHINT1(24)); | |
258 | ||
259 | // TLorentzVector* jet2 = new TLorentzVector(fHijing->GetHINT1(31), | |
260 | // fHijing->GetHINT1(32), | |
261 | // fHijing->GetHINT1(33), | |
262 | // fHijing->GetHINT1(34)); | |
263 | // // After final state gluon radiation. | |
264 | // TLorentzVector* jet3 = new TLorentzVector(fHijing->GetHINT1(26), | |
265 | // fHijing->GetHINT1(27), | |
266 | // fHijing->GetHINT1(28), | |
267 | // fHijing->GetHINT1(29)); | |
268 | ||
269 | // TLorentzVector* jet4 = new TLorentzVector(fHijing->GetHINT1(36), | |
270 | // fHijing->GetHINT1(37), | |
271 | // fHijing->GetHINT1(38), | |
272 | // fHijing->GetHINT1(39)); | |
273 | // ((AliGenHijingEventHeader*) header)->SetJets(jet1, jet2, jet3, jet4); | |
274 | // Bookkeeping for kinematic bias | |
275 | // ((AliGenHijingEventHeader*) header)->SetTrials(fTrials); | |
276 | // Event Vertex | |
277 | header->SetPrimaryVertex(fVertex); | |
278 | AddHeader(header); | |
279 | fCollisionGeometry = (AliGenHijingEventHeader*) header; | |
280 | ||
deb008f3 | 281 | delete idsOnStack; |
282 | ||
c9a371fc | 283 | // gAlice->SetGenEventHeader(header); |
2e967919 | 284 | } |
285 | ||
286 | void AliGenTherminator::Init() | |
287 | { | |
288 | // Initialize global variables and | |
289 | // particle and decay tables | |
290 | if (fFileName.Length() == 0) | |
291 | fFileName = "event.out"; | |
292 | ReadShareParticleTable(); | |
293 | ||
294 | SetMC(new TTherminator()); | |
295 | ||
296 | AliGenMC::Init(); | |
297 | ((TTherminator *) fMCEvGen)->Initialize(); | |
298 | } | |
299 | void AliGenTherminator::SetFileName(const char *infilename) | |
300 | { | |
301 | // Set parameter filename | |
302 | fFileName = infilename; | |
303 | } | |
304 | void AliGenTherminator::SetEventNumberInFile(int evnum) | |
305 | { | |
306 | // Set number of events to generate - default: 1 | |
307 | fEventNumber = evnum; | |
308 | } | |
309 | ||
310 | void AliGenTherminator::ReadShareParticleTable() | |
311 | { | |
312 | // Read in particle table from share | |
313 | // and add missing particle type to TDatabasePDG | |
314 | ||
315 | char str[50]; | |
316 | char str1[200]; | |
317 | ||
318 | TDatabasePDG *tInstance = TDatabasePDG::Instance(); | |
319 | TParticlePDG *tParticleType; | |
320 | ||
321 | AliWarning(Form("Reading particle types from particles.data")); | |
ae89e57a | 322 | |
323 | TString aroot = gSystem->Getenv("ALICE_ROOT"); | |
324 | ifstream in((aroot+"/TTherminator/data/SHARE/particles.data").Data()); | |
325 | // ifstream in("particles.data"); | |
2e967919 | 326 | |
327 | int charge; | |
328 | ||
329 | int number=0; | |
330 | if ((in) && (in.is_open())) | |
331 | { | |
332 | //START OF HEAD-LINE | |
333 | in.ignore(200,'\n'); | |
334 | in.ignore(200,'\n'); | |
335 | in.ignore(200,'\n'); | |
336 | //END OF HEAD-LINE | |
337 | ||
338 | while (in>>str) | |
339 | { | |
340 | if (/*(*str == '#')||*/(*str<65)||(*str>122)) | |
341 | { | |
342 | in.getline(str1,200); | |
343 | continue; | |
344 | } | |
345 | double mass, gamma, spin, tI3, tI, q, s, aq, as, c, ac, mc; | |
346 | ||
347 | in>>mass>>gamma>>spin>>tI>>tI3>>q>>s>>aq>>as>>c>>ac>>mc; | |
348 | number++; | |
349 | tParticleType = tInstance->GetParticle((int) mc); | |
350 | if (!tParticleType) { | |
351 | charge = 0; | |
352 | if (strstr(str, "plu")) charge = 1; | |
353 | if (strstr(str, "min")) charge = -1; | |
354 | if (strstr(str, "plb")) charge = -1; | |
355 | if (strstr(str, "mnb")) charge = 1; | |
356 | if (strstr(str, "plp")) charge = 2; | |
357 | if (strstr(str, "ppb")) charge = -2; | |
358 | tInstance->AddParticle(str, str, mass, gamma == 0.0 ? 1:0, gamma, charge , "meson", (int) mc); | |
359 | AliWarning(Form("Added particle %s with PDG PID %d charge %d", str, (int) mc, charge)); | |
360 | // 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)); | |
361 | ||
362 | } | |
363 | } | |
364 | in.close(); | |
365 | } | |
366 | CreateTherminatorInputFile(); | |
367 | } | |
368 | ||
369 | void AliGenTherminator::CreateTherminatorInputFile() | |
370 | { | |
371 | // Create Therminator input file | |
ae89e57a | 372 | const char *aroot = gSystem->Getenv("ALICE_ROOT"); |
2e967919 | 373 | ofstream *ostr = new ofstream("therminator.in"); |
374 | (*ostr) << "NumberOfEvents = 1" << endl; | |
375 | (*ostr) << "Randomize = 1" << endl; | |
376 | (*ostr) << "TableType = SHARE" << endl; | |
ae89e57a | 377 | (*ostr) << "InputDirSHARE = "<< aroot << "/TTherminator/data/SHARE" << endl; |
2e967919 | 378 | (*ostr) << "EventOutputFile = " << fFileName.Data() << endl; |
ae89e57a | 379 | (*ostr) << "FOHSLocation = " << fFOHSlocation.Data() << endl; |
2e967919 | 380 | (*ostr) << "FreezeOutModel = " << fFreezeOutModel.Data() << endl; |
381 | (*ostr) << "BWVt = " << fBWVt << endl; | |
382 | (*ostr) << "Tau = " << fTau << endl; | |
383 | (*ostr) << "RhoMax = " << fRhoMax << endl; | |
384 | (*ostr) << "Temperature = " << fTemperature << endl; | |
385 | (*ostr) << "MiuI = " << fMiuI << endl; | |
386 | (*ostr) << "MiuS = " << fMiuS << endl; | |
387 | (*ostr) << "MiuB = " << fMiuB << endl; | |
388 | (*ostr) << "AlphaRange = " << fAlfaRange << endl; | |
389 | (*ostr) << "RapidityRange = " << fRapRange << endl; | |
390 | (*ostr) << "NumberOfIntegrateSamples = 1000000" << endl; | |
391 | } | |
392 | ||
393 | void AliGenTherminator::SetModel(const char *model) | |
394 | { | |
395 | // Set the freeze-out model to use | |
396 | fFreezeOutModel = model; | |
397 | } | |
ae89e57a | 398 | |
399 | void AliGenTherminator::SetLhyquidSet(const char *set) | |
400 | { | |
401 | // Select one of pregenerated Lhyquid hypersurfaces | |
402 | const char *aroot = gSystem->Getenv("ALICE_ROOT"); | |
403 | if (strstr(set, "LHC500C0005")) { | |
404 | AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface")); | |
7d4c52b2 | 405 | AliWarning(Form(" Pb-Pb collisions, centrality 0-5 percent")); |
ae89e57a | 406 | AliWarning(Form(" initial temperature at tau=1 fm in the center Ti=500 MeV")); |
407 | AliWarning(Form(" freeze-out criteria Tf=145 MeV")); | |
408 | AliWarning(Form(" for details see $(ALICE_ROOT)/TTherminator/data/LHC500C0005/FO.txt")); | |
409 | fFOHSlocation = aroot; | |
410 | fFOHSlocation += "/TTherminator/data/LHC500C0005"; | |
411 | } | |
412 | else if (strstr(set, "LHC500C2030")) { | |
413 | AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface")); | |
7d4c52b2 | 414 | AliWarning(Form(" Pb-Pb collisions, centrality 20-30 percent")); |
ae89e57a | 415 | AliWarning(Form(" initial temperature at tau=1 fm in the center Ti=500 MeV")); |
416 | AliWarning(Form(" freeze-out criteria Tf=145 MeV")); | |
417 | AliWarning(Form(" for details see $(ALICE_ROOT)/TTherminator/data/LHC500C2030/FO.txt")); | |
418 | fFOHSlocation = aroot; | |
419 | fFOHSlocation += "/TTherminator/data/LHC500C2030"; | |
420 | } | |
421 | else { | |
422 | AliWarning(Form("Did not find Lhyquid set %s", set)); | |
423 | AliWarning(Form("Reverting to default: current directory")); | |
424 | fFOHSlocation = ""; | |
425 | } | |
426 | } | |
427 | ||
428 | void AliGenTherminator::SetLhyquidInputDir(const char *inputdir) | |
429 | { | |
430 | // Select Your own Lhyquid hypersurface | |
431 | fFOHSlocation = inputdir; | |
432 | } |