]>
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 | |
47 | fParticles = new TClonesArray("TParticle",20000); | |
48 | } | |
49 | AliGenTherminator::AliGenTherminator(Int_t npart): | |
50 | AliGenMC(npart), | |
51 | fNt(0), | |
52 | fEventNumber(0), | |
53 | fFileName(""), | |
54 | fFreezeOutModel(""), | |
55 | fTemperature(0.1656), | |
56 | fMiuI(-0.0009), | |
57 | fMiuS(0.0), | |
58 | fMiuB(0.0008), | |
59 | fAlfaRange(6.0), | |
60 | fRapRange(4.0), | |
61 | fRhoMax(8.0), | |
62 | fTau(8.0), | |
63 | fBWA(0.0), | |
64 | fBWVt(1.41), | |
65 | fBWDelay(0.0) | |
66 | { | |
67 | // Constructor specifying the size of the particle table | |
68 | fParticles = new TClonesArray("TParticle",20000); | |
69 | fNprimaries = 0; | |
70 | } | |
71 | ||
72 | AliGenTherminator::~AliGenTherminator() | |
73 | { | |
74 | // AliGenMC::~AliGenMC(); | |
75 | // if (fTherminator) delete fTherminator; | |
76 | } | |
77 | ||
78 | void AliGenTherminator::Generate() | |
79 | { | |
80 | // Run single event generation with the Therminator model | |
81 | AliWarning("Generating event from AliGenTherminator"); | |
82 | ||
83 | Float_t polar[3] = {0,0,0}; | |
84 | Float_t origin[3] = {0,0,0}; | |
85 | Float_t origin0[3] = {0,0,0}; | |
86 | Float_t p[3]; | |
87 | Float_t mass, energy; | |
88 | ||
89 | Int_t nt = 0; | |
90 | Int_t j, kf, ks, imo; | |
91 | kf = 0; | |
92 | ||
93 | Vertex(); | |
94 | for (j=0; j < 3; j++) origin0[j] = fVertex[j]; | |
95 | ||
96 | // Generate one event | |
97 | ||
98 | ((TTherminator *) fMCEvGen)->GenerateEvent(); | |
99 | AliWarning("Generated"); | |
100 | ((TTherminator *) fMCEvGen)->ImportParticles(fParticles); | |
101 | ||
102 | Int_t np = fParticles->GetEntriesFast(); | |
103 | AliWarning(Form("Imported %d particles", np)); | |
104 | ||
105 | TParticle *iparticle; | |
9abb118d | 106 | Double_t evrot = gRandom->Rndm()*TMath::Pi(); |
2e967919 | 107 | |
108 | for (int i = 0; i < np; i++) { | |
109 | iparticle = (TParticle *) fParticles->At(i); | |
110 | Bool_t hasMother = (iparticle->GetFirstMother() >=0); | |
111 | Bool_t hasDaughter = (iparticle->GetFirstDaughter() >=0); | |
112 | ||
113 | kf = iparticle->GetPdgCode(); | |
114 | ks = iparticle->GetStatusCode(); | |
9abb118d | 115 | Double_t aphi = TMath::ATan2(iparticle->Py(), iparticle->Px()); |
116 | Double_t arho = TMath::Hypot(iparticle->Px(), iparticle->Py()); | |
117 | p[0] = arho*TMath::Cos(aphi + evrot); | |
118 | p[1] = arho*TMath::Sin(aphi + evrot); | |
119 | // p[0] = iparticle->Px(); | |
120 | // p[1] = iparticle->Py(); | |
2e967919 | 121 | p[2] = iparticle->Pz(); |
122 | mass = TDatabasePDG::Instance()->GetParticle(kf)->Mass(); | |
123 | energy = sqrt(mass*mass + p[0]*p[0] + p[1]*p[1] + p[2]*p[2]); | |
124 | origin[0] = origin0[0]+iparticle->Vx(); | |
125 | origin[1] = origin0[1]+iparticle->Vy(); | |
126 | origin[2] = origin0[2]+iparticle->Vz(); | |
127 | ||
128 | imo = -1; | |
129 | TParticle* mother = 0; | |
130 | if (hasMother) { | |
131 | imo = iparticle->GetFirstMother(); | |
132 | mother = (TParticle *) fParticles->At(imo); | |
133 | } // if has mother | |
134 | Bool_t tFlag = (!hasDaughter); | |
135 | ||
136 | printf("Pushing Track %d with status %d mother %d\n", kf, tFlag, imo); | |
137 | PushTrack(tFlag,imo,kf, | |
138 | p[0],p[1],p[2],energy, | |
139 | origin[0],origin[1],origin[2],iparticle->T()*0.197327*1e-13/300000000, | |
140 | polar[0],polar[1],polar[2], | |
141 | hasMother ? kPDecay:kPNoProcess,nt); | |
142 | ||
143 | fNprimaries++; | |
144 | KeepTrack(nt); | |
145 | } | |
146 | ||
147 | SetHighWaterMark(fNprimaries); | |
148 | ||
149 | TArrayF eventVertex; | |
150 | eventVertex.Set(3); | |
151 | eventVertex[0] = origin0[0]; | |
152 | eventVertex[1] = origin0[1]; | |
153 | eventVertex[2] = origin0[2]; | |
154 | ||
9abb118d | 155 | // Builds the event header, to be called after each event |
156 | AliGenEventHeader* header = new AliGenHijingEventHeader("Therminator"); | |
157 | ||
2e967919 | 158 | // Header |
9abb118d | 159 | // AliGenEventHeader* header = new AliGenEventHeader("Therminator"); |
2e967919 | 160 | // Event Vertex |
9abb118d | 161 | // header->SetPrimaryVertex(eventVertex); |
162 | // header->SetNProduced(fNprimaries); | |
163 | ||
164 | ((AliGenHijingEventHeader*) header)->SetNProduced(fNprimaries); | |
165 | ((AliGenHijingEventHeader*) header)->SetPrimaryVertex(eventVertex); | |
166 | ((AliGenHijingEventHeader*) header)->SetImpactParameter(0.0); | |
167 | ((AliGenHijingEventHeader*) header)->SetTotalEnergy(0.0); | |
168 | ((AliGenHijingEventHeader*) header)->SetHardScatters(0); | |
169 | ((AliGenHijingEventHeader*) header)->SetParticipants(0, 0); | |
170 | ((AliGenHijingEventHeader*) header)->SetCollisions(0, 0, 0, 0); | |
171 | ((AliGenHijingEventHeader*) header)->SetSpectators(0, 0, 0, 0); | |
172 | ((AliGenHijingEventHeader*) header)->SetReactionPlaneAngle(evrot); | |
173 | ||
174 | ||
175 | // 4-momentum vectors of the triggered jets. | |
176 | // | |
177 | // Before final state gluon radiation. | |
178 | // TLorentzVector* jet1 = new TLorentzVector(fHijing->GetHINT1(21), | |
179 | // fHijing->GetHINT1(22), | |
180 | // fHijing->GetHINT1(23), | |
181 | // fHijing->GetHINT1(24)); | |
182 | ||
183 | // TLorentzVector* jet2 = new TLorentzVector(fHijing->GetHINT1(31), | |
184 | // fHijing->GetHINT1(32), | |
185 | // fHijing->GetHINT1(33), | |
186 | // fHijing->GetHINT1(34)); | |
187 | // // After final state gluon radiation. | |
188 | // TLorentzVector* jet3 = new TLorentzVector(fHijing->GetHINT1(26), | |
189 | // fHijing->GetHINT1(27), | |
190 | // fHijing->GetHINT1(28), | |
191 | // fHijing->GetHINT1(29)); | |
192 | ||
193 | // TLorentzVector* jet4 = new TLorentzVector(fHijing->GetHINT1(36), | |
194 | // fHijing->GetHINT1(37), | |
195 | // fHijing->GetHINT1(38), | |
196 | // fHijing->GetHINT1(39)); | |
197 | // ((AliGenHijingEventHeader*) header)->SetJets(jet1, jet2, jet3, jet4); | |
198 | // Bookkeeping for kinematic bias | |
199 | // ((AliGenHijingEventHeader*) header)->SetTrials(fTrials); | |
200 | // Event Vertex | |
201 | header->SetPrimaryVertex(fVertex); | |
202 | AddHeader(header); | |
203 | fCollisionGeometry = (AliGenHijingEventHeader*) header; | |
204 | ||
2e967919 | 205 | gAlice->SetGenEventHeader(header); |
206 | } | |
207 | ||
208 | void AliGenTherminator::Init() | |
209 | { | |
210 | // Initialize global variables and | |
211 | // particle and decay tables | |
212 | if (fFileName.Length() == 0) | |
213 | fFileName = "event.out"; | |
214 | ReadShareParticleTable(); | |
215 | ||
216 | SetMC(new TTherminator()); | |
217 | ||
218 | AliGenMC::Init(); | |
219 | ((TTherminator *) fMCEvGen)->Initialize(); | |
220 | } | |
221 | void AliGenTherminator::SetFileName(const char *infilename) | |
222 | { | |
223 | // Set parameter filename | |
224 | fFileName = infilename; | |
225 | } | |
226 | void AliGenTherminator::SetEventNumberInFile(int evnum) | |
227 | { | |
228 | // Set number of events to generate - default: 1 | |
229 | fEventNumber = evnum; | |
230 | } | |
231 | ||
232 | void AliGenTherminator::ReadShareParticleTable() | |
233 | { | |
234 | // Read in particle table from share | |
235 | // and add missing particle type to TDatabasePDG | |
236 | ||
237 | char str[50]; | |
238 | char str1[200]; | |
239 | ||
240 | TDatabasePDG *tInstance = TDatabasePDG::Instance(); | |
241 | TParticlePDG *tParticleType; | |
242 | ||
243 | AliWarning(Form("Reading particle types from particles.data")); | |
ae89e57a | 244 | |
245 | TString aroot = gSystem->Getenv("ALICE_ROOT"); | |
246 | ifstream in((aroot+"/TTherminator/data/SHARE/particles.data").Data()); | |
247 | // ifstream in("particles.data"); | |
2e967919 | 248 | |
249 | int charge; | |
250 | ||
251 | int number=0; | |
252 | if ((in) && (in.is_open())) | |
253 | { | |
254 | //START OF HEAD-LINE | |
255 | in.ignore(200,'\n'); | |
256 | in.ignore(200,'\n'); | |
257 | in.ignore(200,'\n'); | |
258 | //END OF HEAD-LINE | |
259 | ||
260 | while (in>>str) | |
261 | { | |
262 | if (/*(*str == '#')||*/(*str<65)||(*str>122)) | |
263 | { | |
264 | in.getline(str1,200); | |
265 | continue; | |
266 | } | |
267 | double mass, gamma, spin, tI3, tI, q, s, aq, as, c, ac, mc; | |
268 | ||
269 | in>>mass>>gamma>>spin>>tI>>tI3>>q>>s>>aq>>as>>c>>ac>>mc; | |
270 | number++; | |
271 | tParticleType = tInstance->GetParticle((int) mc); | |
272 | if (!tParticleType) { | |
273 | charge = 0; | |
274 | if (strstr(str, "plu")) charge = 1; | |
275 | if (strstr(str, "min")) charge = -1; | |
276 | if (strstr(str, "plb")) charge = -1; | |
277 | if (strstr(str, "mnb")) charge = 1; | |
278 | if (strstr(str, "plp")) charge = 2; | |
279 | if (strstr(str, "ppb")) charge = -2; | |
280 | tInstance->AddParticle(str, str, mass, gamma == 0.0 ? 1:0, gamma, charge , "meson", (int) mc); | |
281 | AliWarning(Form("Added particle %s with PDG PID %d charge %d", str, (int) mc, charge)); | |
282 | // 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)); | |
283 | ||
284 | } | |
285 | } | |
286 | in.close(); | |
287 | } | |
288 | CreateTherminatorInputFile(); | |
289 | } | |
290 | ||
291 | void AliGenTherminator::CreateTherminatorInputFile() | |
292 | { | |
293 | // Create Therminator input file | |
ae89e57a | 294 | const char *aroot = gSystem->Getenv("ALICE_ROOT"); |
2e967919 | 295 | ofstream *ostr = new ofstream("therminator.in"); |
296 | (*ostr) << "NumberOfEvents = 1" << endl; | |
297 | (*ostr) << "Randomize = 1" << endl; | |
298 | (*ostr) << "TableType = SHARE" << endl; | |
ae89e57a | 299 | (*ostr) << "InputDirSHARE = "<< aroot << "/TTherminator/data/SHARE" << endl; |
2e967919 | 300 | (*ostr) << "EventOutputFile = " << fFileName.Data() << endl; |
ae89e57a | 301 | (*ostr) << "FOHSLocation = " << fFOHSlocation.Data() << endl; |
2e967919 | 302 | (*ostr) << "FreezeOutModel = " << fFreezeOutModel.Data() << endl; |
303 | (*ostr) << "BWVt = " << fBWVt << endl; | |
304 | (*ostr) << "Tau = " << fTau << endl; | |
305 | (*ostr) << "RhoMax = " << fRhoMax << endl; | |
306 | (*ostr) << "Temperature = " << fTemperature << endl; | |
307 | (*ostr) << "MiuI = " << fMiuI << endl; | |
308 | (*ostr) << "MiuS = " << fMiuS << endl; | |
309 | (*ostr) << "MiuB = " << fMiuB << endl; | |
310 | (*ostr) << "AlphaRange = " << fAlfaRange << endl; | |
311 | (*ostr) << "RapidityRange = " << fRapRange << endl; | |
312 | (*ostr) << "NumberOfIntegrateSamples = 1000000" << endl; | |
313 | } | |
314 | ||
315 | void AliGenTherminator::SetModel(const char *model) | |
316 | { | |
317 | // Set the freeze-out model to use | |
318 | fFreezeOutModel = model; | |
319 | } | |
ae89e57a | 320 | |
321 | void AliGenTherminator::SetLhyquidSet(const char *set) | |
322 | { | |
323 | // Select one of pregenerated Lhyquid hypersurfaces | |
324 | const char *aroot = gSystem->Getenv("ALICE_ROOT"); | |
325 | if (strstr(set, "LHC500C0005")) { | |
326 | AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface")); | |
327 | AliWarning(Form(" Pb-Pb collisions, centrality 0-5%")); | |
328 | AliWarning(Form(" initial temperature at tau=1 fm in the center Ti=500 MeV")); | |
329 | AliWarning(Form(" freeze-out criteria Tf=145 MeV")); | |
330 | AliWarning(Form(" for details see $(ALICE_ROOT)/TTherminator/data/LHC500C0005/FO.txt")); | |
331 | fFOHSlocation = aroot; | |
332 | fFOHSlocation += "/TTherminator/data/LHC500C0005"; | |
333 | } | |
334 | else if (strstr(set, "LHC500C2030")) { | |
335 | AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface")); | |
336 | AliWarning(Form(" Pb-Pb collisions, centrality 20-30%")); | |
337 | AliWarning(Form(" initial temperature at tau=1 fm in the center Ti=500 MeV")); | |
338 | AliWarning(Form(" freeze-out criteria Tf=145 MeV")); | |
339 | AliWarning(Form(" for details see $(ALICE_ROOT)/TTherminator/data/LHC500C2030/FO.txt")); | |
340 | fFOHSlocation = aroot; | |
341 | fFOHSlocation += "/TTherminator/data/LHC500C2030"; | |
342 | } | |
343 | else { | |
344 | AliWarning(Form("Did not find Lhyquid set %s", set)); | |
345 | AliWarning(Form("Reverting to default: current directory")); | |
346 | fFOHSlocation = ""; | |
347 | } | |
348 | } | |
349 | ||
350 | void AliGenTherminator::SetLhyquidInputDir(const char *inputdir) | |
351 | { | |
352 | // Select Your own Lhyquid hypersurface | |
353 | fFOHSlocation = inputdir; | |
354 | } |