1 // ALICE event generator based on the THERMINATOR model
2 // It reads the test output of the model and puts it onto
4 // It has an option to use the Lhyquid3D input freeze-out
6 // Author: Adam.Kisiel@cern.ch
11 #include <TClonesArray.h>
12 #include <TMCProcess.h>
13 #include <TDatabasePDG.h>
14 #include <TParticle.h>
17 #include "AliDecayer.h"
18 #include "AliGenEventHeader.h"
19 #include "AliGenHijingEventHeader.h"
20 #include "AliGenTherminator.h"
24 ClassImp(AliGenTherminator)
28 AliGenTherminator::AliGenTherminator():
46 // Default constructor
48 AliGenTherminator::AliGenTherminator(Int_t npart):
66 // Constructor specifying the size of the particle table
70 AliGenTherminator::~AliGenTherminator()
72 // AliGenMC::~AliGenMC();
73 // if (fTherminator) delete fTherminator;
76 void AliGenTherminator::Generate()
78 // Run single event generation with the Therminator model
79 AliWarning("Generating event from AliGenTherminator");
81 Float_t polar[3] = {0,0,0};
82 Float_t origin[3] = {0,0,0};
83 Float_t origin0[3] = {0,0,0};
92 for (j=0; j < 3; j++) origin0[j] = fVertex[j];
96 ((TTherminator *) fMCEvGen)->GenerateEvent();
97 AliWarning("Generated");
98 ((TTherminator *) fMCEvGen)->ImportParticles(&fParticles);
100 Int_t np = fParticles.GetEntriesFast();
101 AliWarning(Form("Imported %d particles", np));
103 TParticle *iparticle;
104 Double_t evrot = gRandom->Rndm()*TMath::Pi();
106 for (int i = 0; i < np; i++) {
107 iparticle = (TParticle *) fParticles.At(i);
108 Bool_t hasMother = (iparticle->GetFirstMother() >=0);
109 Bool_t hasDaughter = (iparticle->GetFirstDaughter() >=0);
112 // This particle has decayed
113 // It will not be tracked
114 // Add it only once with coorduinates not
115 // smeared with primary vertex position
117 kf = iparticle->GetPdgCode();
118 ks = iparticle->GetStatusCode();
119 Double_t aphi = TMath::ATan2(iparticle->Py(), iparticle->Px());
120 Double_t arho = TMath::Hypot(iparticle->Px(), iparticle->Py());
121 p[0] = arho*TMath::Cos(aphi + evrot);
122 p[1] = arho*TMath::Sin(aphi + evrot);
123 // p[0] = iparticle->Px();
124 // p[1] = iparticle->Py();
125 p[2] = iparticle->Pz();
126 mass = TDatabasePDG::Instance()->GetParticle(kf)->Mass();
127 energy = sqrt(mass*mass + p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
129 Double_t vphi = TMath::ATan2(iparticle->Vy(), iparticle->Vx());
130 Double_t vrho = TMath::Hypot(iparticle->Vx(), iparticle->Vy());
131 origin[0] = origin0[0]+vrho*TMath::Cos(vphi + evrot);
132 origin[1] = origin0[1]+vrho*TMath::Sin(vphi + evrot);
133 origin[2] = origin0[2]+iparticle->Vz();
136 TParticle* mother = 0;
138 imo = iparticle->GetFirstMother();
139 mother = (TParticle *) fParticles->At(imo);
141 Bool_t tFlag = (!hasDaughter);
143 printf("Pushing Track %d with status %d mother %d\n", kf, tFlag, imo);
144 PushTrack(tFlag,imo,kf,
145 p[0],p[1],p[2],energy,
146 origin[0],origin[1],origin[2],iparticle->T()*0.197327*1e-13/300000000,
147 polar[0],polar[1],polar[2],
148 hasMother ? kPDecay:kPNoProcess,nt);
153 // This is a final state particle
154 // It will be tracked
155 // Add it TWICE to the stack !!!
156 // First time with event-wide coordicates (for femtoscopy) -
157 // this one will not be tracked
158 // Second time with event-wide ccordiantes and vertex smearing
159 // this one will be tracked
161 kf = iparticle->GetPdgCode();
162 ks = iparticle->GetStatusCode();
163 Double_t aphi = TMath::ATan2(iparticle->Py(), iparticle->Px());
164 Double_t arho = TMath::Hypot(iparticle->Px(), iparticle->Py());
165 p[0] = arho*TMath::Cos(aphi + evrot);
166 p[1] = arho*TMath::Sin(aphi + evrot);
167 // p[0] = iparticle->Px();
168 // p[1] = iparticle->Py();
169 p[2] = iparticle->Pz();
170 mass = TDatabasePDG::Instance()->GetParticle(kf)->Mass();
171 energy = sqrt(mass*mass + p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
173 Double_t vphi = TMath::ATan2(iparticle->Vy(), iparticle->Vx());
174 Double_t vrho = TMath::Hypot(iparticle->Vx(), iparticle->Vy());
175 origin[0] = vrho*TMath::Cos(vphi + evrot);
176 origin[1] = vrho*TMath::Sin(vphi + evrot);
177 origin[2] = iparticle->Vz();
180 TParticle* mother = 0;
182 imo = iparticle->GetFirstMother();
183 mother = (TParticle *) fParticles->At(imo);
185 Bool_t tFlag = (hasDaughter);
187 printf("Pushing Track %d with status %d mother %d\n", kf, tFlag, imo);
188 PushTrack(tFlag,imo,kf,
189 p[0],p[1],p[2],energy,
190 origin[0],origin[1],origin[2],iparticle->T()*0.197327*1e-13/300000000,
191 polar[0],polar[1],polar[2],
192 hasMother ? kPDecay:kPNoProcess,nt);
196 origin[0] = origin0[0]+vrho*TMath::Cos(vphi + evrot);
197 origin[1] = origin0[1]+vrho*TMath::Sin(vphi + evrot);
198 origin[2] = origin0[2]+iparticle->Vz();
201 mother = (TParticle *) fParticles->At(nt);
202 tFlag = (!hasDaughter);
204 printf("Pushing Track %d with status %d mother %d\n", kf, tFlag, imo);
205 PushTrack(tFlag,imo,kf,
206 p[0],p[1],p[2],energy,
207 origin[0],origin[1],origin[2],iparticle->T()*0.197327*1e-13/300000000,
208 polar[0],polar[1],polar[2],
209 hasMother ? kPDecay:kPNoProcess,nt);
215 SetHighWaterMark(fNprimaries);
219 eventVertex[0] = origin0[0];
220 eventVertex[1] = origin0[1];
221 eventVertex[2] = origin0[2];
223 // Builds the event header, to be called after each event
224 AliGenEventHeader* header = new AliGenHijingEventHeader("Therminator");
227 // AliGenEventHeader* header = new AliGenEventHeader("Therminator");
229 // header->SetPrimaryVertex(eventVertex);
230 // header->SetNProduced(fNprimaries);
232 ((AliGenHijingEventHeader*) header)->SetNProduced(fNprimaries);
233 ((AliGenHijingEventHeader*) header)->SetPrimaryVertex(eventVertex);
234 ((AliGenHijingEventHeader*) header)->SetImpactParameter(0.0);
235 ((AliGenHijingEventHeader*) header)->SetTotalEnergy(0.0);
236 ((AliGenHijingEventHeader*) header)->SetHardScatters(0);
237 ((AliGenHijingEventHeader*) header)->SetParticipants(0, 0);
238 ((AliGenHijingEventHeader*) header)->SetCollisions(0, 0, 0, 0);
239 ((AliGenHijingEventHeader*) header)->SetSpectators(0, 0, 0, 0);
240 ((AliGenHijingEventHeader*) header)->SetReactionPlaneAngle(evrot);
243 // 4-momentum vectors of the triggered jets.
245 // Before final state gluon radiation.
246 // TLorentzVector* jet1 = new TLorentzVector(fHijing->GetHINT1(21),
247 // fHijing->GetHINT1(22),
248 // fHijing->GetHINT1(23),
249 // fHijing->GetHINT1(24));
251 // TLorentzVector* jet2 = new TLorentzVector(fHijing->GetHINT1(31),
252 // fHijing->GetHINT1(32),
253 // fHijing->GetHINT1(33),
254 // fHijing->GetHINT1(34));
255 // // After final state gluon radiation.
256 // TLorentzVector* jet3 = new TLorentzVector(fHijing->GetHINT1(26),
257 // fHijing->GetHINT1(27),
258 // fHijing->GetHINT1(28),
259 // fHijing->GetHINT1(29));
261 // TLorentzVector* jet4 = new TLorentzVector(fHijing->GetHINT1(36),
262 // fHijing->GetHINT1(37),
263 // fHijing->GetHINT1(38),
264 // fHijing->GetHINT1(39));
265 // ((AliGenHijingEventHeader*) header)->SetJets(jet1, jet2, jet3, jet4);
266 // Bookkeeping for kinematic bias
267 // ((AliGenHijingEventHeader*) header)->SetTrials(fTrials);
269 header->SetPrimaryVertex(fVertex);
271 fCollisionGeometry = (AliGenHijingEventHeader*) header;
273 // gAlice->SetGenEventHeader(header);
276 void AliGenTherminator::Init()
278 // Initialize global variables and
279 // particle and decay tables
280 if (fFileName.Length() == 0)
281 fFileName = "event.out";
282 ReadShareParticleTable();
284 SetMC(new TTherminator());
287 ((TTherminator *) fMCEvGen)->Initialize();
289 void AliGenTherminator::SetFileName(const char *infilename)
291 // Set parameter filename
292 fFileName = infilename;
294 void AliGenTherminator::SetEventNumberInFile(int evnum)
296 // Set number of events to generate - default: 1
297 fEventNumber = evnum;
300 void AliGenTherminator::ReadShareParticleTable()
302 // Read in particle table from share
303 // and add missing particle type to TDatabasePDG
308 TDatabasePDG *tInstance = TDatabasePDG::Instance();
309 TParticlePDG *tParticleType;
311 AliWarning(Form("Reading particle types from particles.data"));
313 TString aroot = gSystem->Getenv("ALICE_ROOT");
314 ifstream in((aroot+"/TTherminator/data/SHARE/particles.data").Data());
315 // ifstream in("particles.data");
320 if ((in) && (in.is_open()))
330 if (/*(*str == '#')||*/(*str<65)||(*str>122))
332 in.getline(str1,200);
335 double mass, gamma, spin, tI3, tI, q, s, aq, as, c, ac, mc;
337 in>>mass>>gamma>>spin>>tI>>tI3>>q>>s>>aq>>as>>c>>ac>>mc;
339 tParticleType = tInstance->GetParticle((int) mc);
340 if (!tParticleType) {
342 if (strstr(str, "plu")) charge = 1;
343 if (strstr(str, "min")) charge = -1;
344 if (strstr(str, "plb")) charge = -1;
345 if (strstr(str, "mnb")) charge = 1;
346 if (strstr(str, "plp")) charge = 2;
347 if (strstr(str, "ppb")) charge = -2;
348 tInstance->AddParticle(str, str, mass, gamma == 0.0 ? 1:0, gamma, charge , "meson", (int) mc);
349 AliWarning(Form("Added particle %s with PDG PID %d charge %d", str, (int) mc, charge));
350 // 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));
356 CreateTherminatorInputFile();
359 void AliGenTherminator::CreateTherminatorInputFile()
361 // Create Therminator input file
362 const char *aroot = gSystem->Getenv("ALICE_ROOT");
363 ofstream *ostr = new ofstream("therminator.in");
364 (*ostr) << "NumberOfEvents = 1" << endl;
365 (*ostr) << "Randomize = 1" << endl;
366 (*ostr) << "TableType = SHARE" << endl;
367 (*ostr) << "InputDirSHARE = "<< aroot << "/TTherminator/data/SHARE" << endl;
368 (*ostr) << "EventOutputFile = " << fFileName.Data() << endl;
369 (*ostr) << "FOHSLocation = " << fFOHSlocation.Data() << endl;
370 (*ostr) << "FreezeOutModel = " << fFreezeOutModel.Data() << endl;
371 (*ostr) << "BWVt = " << fBWVt << endl;
372 (*ostr) << "Tau = " << fTau << endl;
373 (*ostr) << "RhoMax = " << fRhoMax << endl;
374 (*ostr) << "Temperature = " << fTemperature << endl;
375 (*ostr) << "MiuI = " << fMiuI << endl;
376 (*ostr) << "MiuS = " << fMiuS << endl;
377 (*ostr) << "MiuB = " << fMiuB << endl;
378 (*ostr) << "AlphaRange = " << fAlfaRange << endl;
379 (*ostr) << "RapidityRange = " << fRapRange << endl;
380 (*ostr) << "NumberOfIntegrateSamples = 1000000" << endl;
383 void AliGenTherminator::SetModel(const char *model)
385 // Set the freeze-out model to use
386 fFreezeOutModel = model;
389 void AliGenTherminator::SetLhyquidSet(const char *set)
391 // Select one of pregenerated Lhyquid hypersurfaces
392 const char *aroot = gSystem->Getenv("ALICE_ROOT");
393 if (strstr(set, "LHC500C0005")) {
394 AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
395 AliWarning(Form(" Pb-Pb collisions, centrality 0-5%"));
396 AliWarning(Form(" initial temperature at tau=1 fm in the center Ti=500 MeV"));
397 AliWarning(Form(" freeze-out criteria Tf=145 MeV"));
398 AliWarning(Form(" for details see $(ALICE_ROOT)/TTherminator/data/LHC500C0005/FO.txt"));
399 fFOHSlocation = aroot;
400 fFOHSlocation += "/TTherminator/data/LHC500C0005";
402 else if (strstr(set, "LHC500C2030")) {
403 AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
404 AliWarning(Form(" Pb-Pb collisions, centrality 20-30%"));
405 AliWarning(Form(" initial temperature at tau=1 fm in the center Ti=500 MeV"));
406 AliWarning(Form(" freeze-out criteria Tf=145 MeV"));
407 AliWarning(Form(" for details see $(ALICE_ROOT)/TTherminator/data/LHC500C2030/FO.txt"));
408 fFOHSlocation = aroot;
409 fFOHSlocation += "/TTherminator/data/LHC500C2030";
412 AliWarning(Form("Did not find Lhyquid set %s", set));
413 AliWarning(Form("Reverting to default: current directory"));
418 void AliGenTherminator::SetLhyquidInputDir(const char *inputdir)
420 // Select Your own Lhyquid hypersurface
421 fFOHSlocation = inputdir;