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():
47 // Default constructor
58 AliGenTherminator::AliGenTherminator(Int_t npart):
77 // Constructor specifying the size of the particle table
88 AliGenTherminator::~AliGenTherminator()
90 // AliGenMC::~AliGenMC();
91 // if (fTherminator) delete fTherminator;
94 void AliGenTherminator::Generate()
96 // Run single event generation with the Therminator model
97 AliWarning("Generating event from AliGenTherminator");
99 Float_t polar[3] = {0,0,0};
100 Float_t origin[3] = {0,0,0};
101 Float_t origin0[3] = {0,0,0};
103 Float_t mass, energy;
106 Int_t j, kf, ks, imo;
110 for (j=0; j < 3; j++) origin0[j] = fVertex[j];
112 // Generate one event
114 ((TTherminator *) fMCEvGen)->GenerateEvent();
115 AliWarning("Generated");
116 ((TTherminator *) fMCEvGen)->ImportParticles(&fParticles);
118 Int_t np = fParticles.GetEntriesFast();
119 AliWarning(Form("Imported %d particles", np));
122 idsOnStack = new Int_t[np];
124 TParticle *iparticle;
125 Double_t evrot = gRandom->Rndm()*TMath::Pi();
127 for (int i = 0; i < np; i++) {
128 iparticle = (TParticle *) fParticles.At(i);
129 Bool_t hasMother = (iparticle->GetFirstMother() >=0);
130 Bool_t hasDaughter = (iparticle->GetFirstDaughter() >=0);
133 // This particle has decayed
134 // It will not be tracked
135 // Add it only once with coorduinates not
136 // smeared with primary vertex position
138 kf = iparticle->GetPdgCode();
139 ks = iparticle->GetStatusCode();
140 Double_t aphi = TMath::ATan2(iparticle->Py(), iparticle->Px());
141 Double_t arho = TMath::Hypot(iparticle->Px(), iparticle->Py());
142 p[0] = arho*TMath::Cos(aphi + evrot);
143 p[1] = arho*TMath::Sin(aphi + evrot);
144 // p[0] = iparticle->Px();
145 // p[1] = iparticle->Py();
146 p[2] = iparticle->Pz();
147 mass = TDatabasePDG::Instance()->GetParticle(kf)->Mass();
148 energy = sqrt(mass*mass + p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
150 Double_t vphi = TMath::ATan2(iparticle->Vy(), iparticle->Vx());
151 Double_t vrho = TMath::Hypot(iparticle->Vx(), iparticle->Vy());
152 origin[0] = vrho*TMath::Cos(vphi + evrot);
153 origin[1] = vrho*TMath::Sin(vphi + evrot);
154 origin[2] = iparticle->Vz();
157 // TParticle* mother = 0;
159 imo = iparticle->GetFirstMother();
160 // mother = (TParticle *) fParticles.At(imo);
162 Bool_t tFlag = (!hasDaughter);
164 // printf("Pushing Track %d with status %d mother %d\n", kf, tFlag, imo>=0?idsOnStack[imo]:imo);
165 PushTrack(tFlag,imo>=0?idsOnStack[imo]:imo,kf,
166 p[0],p[1],p[2],energy,
167 origin[0],origin[1],origin[2],iparticle->T(),
168 polar[0],polar[1],polar[2],
169 hasMother ? kPDecay:kPNoProcess,nt);
175 // This is a final state particle
176 // It will be tracked
177 // Add it TWICE to the stack !!!
178 // First time with event-wide coordicates (for femtoscopy) -
179 // this one will not be tracked
180 // Second time with event-wide ccordiantes and vertex smearing
181 // this one will be tracked
183 kf = iparticle->GetPdgCode();
184 ks = iparticle->GetStatusCode();
185 Double_t aphi = TMath::ATan2(iparticle->Py(), iparticle->Px());
186 Double_t arho = TMath::Hypot(iparticle->Px(), iparticle->Py());
187 p[0] = arho*TMath::Cos(aphi + evrot);
188 p[1] = arho*TMath::Sin(aphi + evrot);
189 // p[0] = iparticle->Px();
190 // p[1] = iparticle->Py();
191 p[2] = iparticle->Pz();
192 mass = TDatabasePDG::Instance()->GetParticle(kf)->Mass();
193 energy = sqrt(mass*mass + p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
195 Double_t vphi = TMath::ATan2(iparticle->Vy(), iparticle->Vx());
196 Double_t vrho = TMath::Hypot(iparticle->Vx(), iparticle->Vy());
197 origin[0] = vrho*TMath::Cos(vphi + evrot);
198 origin[1] = vrho*TMath::Sin(vphi + evrot);
199 origin[2] = iparticle->Vz();
202 // TParticle* mother = 0;
204 imo = iparticle->GetFirstMother();
205 // mother = (TParticle *) fParticles.At(imo);
207 Bool_t tFlag = (hasDaughter);
209 // printf("Found mother %i with true id %i\n", imo, imo>=0?idsOnStack[imo]:imo);
210 // printf("Pushing Track %d with status %d mother %d\n", kf, tFlag, imo>=0?idsOnStack[imo]:imo);
211 PushTrack(tFlag,imo>=0?idsOnStack[imo]:imo,kf,
212 p[0],p[1],p[2],energy,
213 origin[0],origin[1],origin[2],iparticle->T(),
214 polar[0],polar[1],polar[2],
215 hasMother ? kPDecay:kPNoProcess,nt);
220 origin[0] = origin0[0]+vrho*TMath::Cos(vphi + evrot);
221 origin[1] = origin0[1]+vrho*TMath::Sin(vphi + evrot);
222 origin[2] = origin0[2]+iparticle->Vz();
225 // mother = (TParticle *) fParticles.At(nt);
226 tFlag = (!hasDaughter);
228 // printf("Pushing Track %d with status %d mother %d\n", kf, tFlag, imo);
229 PushTrack(tFlag,imo,kf,
230 p[0],p[1],p[2],energy,
231 origin[0],origin[1],origin[2],iparticle->T(),
232 polar[0],polar[1],polar[2],
233 hasMother ? kPDecay:kPNoProcess,nt);
241 SetHighWaterMark(fNprimaries);
245 eventVertex[0] = origin0[0];
246 eventVertex[1] = origin0[1];
247 eventVertex[2] = origin0[2];
249 // Builds the event header, to be called after each event
250 AliGenEventHeader* header = new AliGenHijingEventHeader("Therminator");
253 // AliGenEventHeader* header = new AliGenEventHeader("Therminator");
255 // header->SetPrimaryVertex(eventVertex);
256 // header->SetNProduced(fNprimaries);
258 ((AliGenHijingEventHeader*) header)->SetNProduced(fNprimaries);
259 ((AliGenHijingEventHeader*) header)->SetPrimaryVertex(eventVertex);
260 ((AliGenHijingEventHeader*) header)->SetImpactParameter(0.0);
261 ((AliGenHijingEventHeader*) header)->SetTotalEnergy(0.0);
262 ((AliGenHijingEventHeader*) header)->SetHardScatters(0);
263 ((AliGenHijingEventHeader*) header)->SetParticipants(0, 0);
264 ((AliGenHijingEventHeader*) header)->SetCollisions(0, 0, 0, 0);
265 ((AliGenHijingEventHeader*) header)->SetSpectators(0, 0, 0, 0);
266 ((AliGenHijingEventHeader*) header)->SetReactionPlaneAngle(evrot);
269 // 4-momentum vectors of the triggered jets.
271 // Before final state gluon radiation.
272 // TLorentzVector* jet1 = new TLorentzVector(fHijing->GetHINT1(21),
273 // fHijing->GetHINT1(22),
274 // fHijing->GetHINT1(23),
275 // fHijing->GetHINT1(24));
277 // TLorentzVector* jet2 = new TLorentzVector(fHijing->GetHINT1(31),
278 // fHijing->GetHINT1(32),
279 // fHijing->GetHINT1(33),
280 // fHijing->GetHINT1(34));
281 // // After final state gluon radiation.
282 // TLorentzVector* jet3 = new TLorentzVector(fHijing->GetHINT1(26),
283 // fHijing->GetHINT1(27),
284 // fHijing->GetHINT1(28),
285 // fHijing->GetHINT1(29));
287 // TLorentzVector* jet4 = new TLorentzVector(fHijing->GetHINT1(36),
288 // fHijing->GetHINT1(37),
289 // fHijing->GetHINT1(38),
290 // fHijing->GetHINT1(39));
291 // ((AliGenHijingEventHeader*) header)->SetJets(jet1, jet2, jet3, jet4);
292 // Bookkeeping for kinematic bias
293 // ((AliGenHijingEventHeader*) header)->SetTrials(fTrials);
295 header->SetPrimaryVertex(fVertex);
297 fCollisionGeometry = (AliGenHijingEventHeader*) header;
299 delete [] idsOnStack;
301 // gAlice->SetGenEventHeader(header);
304 void AliGenTherminator::Init()
306 // Initialize global variables and
307 // particle and decay tables
308 if (fFileName.Length() == 0)
309 fFileName = "event.out";
310 ReadShareParticleTable();
312 SetMC(new TTherminator());
315 ((TTherminator *) fMCEvGen)->Initialize();
317 void AliGenTherminator::SetFileName(const char *infilename)
319 // Set parameter filename
320 fFileName = infilename;
322 void AliGenTherminator::SetEventNumberInFile(int evnum)
324 // Set number of events to generate - default: 1
325 fEventNumber = evnum;
328 void AliGenTherminator::ReadShareParticleTable()
330 // Read in particle table from share
331 // and add missing particle type to TDatabasePDG
336 TDatabasePDG *tInstance = TDatabasePDG::Instance();
337 TParticlePDG *tParticleType;
339 AliWarning(Form("Reading particle types from particles.data"));
341 TString aroot = gSystem->Getenv("ALICE_ROOT");
342 ifstream in((aroot+"/TTherminator/data/SHARE/particles.data").Data());
343 // ifstream in("particles.data");
348 if ((in) && (in.is_open()))
358 if (/*(*str == '#')||*/(*str<65)||(*str>122))
360 in.getline(str1,200);
363 double mass, gamma, spin, tI3, tI, q, s, aq, as, c, ac, mc;
365 in>>mass>>gamma>>spin>>tI>>tI3>>q>>s>>aq>>as>>c>>ac>>mc;
367 tParticleType = tInstance->GetParticle((int) mc);
368 if (!tParticleType) {
370 if (strstr(str, "plu")) charge = 1;
371 if (strstr(str, "min")) charge = -1;
372 if (strstr(str, "plb")) charge = -1;
373 if (strstr(str, "mnb")) charge = 1;
374 if (strstr(str, "plp")) charge = 2;
375 if (strstr(str, "ppb")) charge = -2;
376 tInstance->AddParticle(str, str, mass, gamma == 0.0 ? 1:0, gamma, charge , "meson", (int) mc);
377 AliWarning(Form("Added particle %s with PDG PID %d charge %d", str, (int) mc, charge));
378 // 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));
384 CreateTherminatorInputFile();
387 void AliGenTherminator::CreateTherminatorInputFile()
389 // Create Therminator input file
390 const char *aroot = gSystem->Getenv("ALICE_ROOT");
391 ofstream *ostr = new ofstream("therminator.in");
392 (*ostr) << "NumberOfEvents = 1" << endl;
393 (*ostr) << "Randomize = 1" << endl;
394 (*ostr) << "TableType = SHARE" << endl;
395 (*ostr) << "InputDirSHARE = "<< aroot << "/TTherminator/data/SHARE" << endl;
396 (*ostr) << "EventOutputFile = " << fFileName.Data() << endl;
397 (*ostr) << "FOHSLocation = " << fFOHSlocation.Data() << endl;
398 (*ostr) << "FreezeOutModel = " << fFreezeOutModel.Data() << endl;
399 (*ostr) << "BWVt = " << fBWVt << endl;
400 (*ostr) << "Tau = " << fTau << endl;
401 (*ostr) << "RhoMax = " << fRhoMax << endl;
402 (*ostr) << "Temperature = " << fTemperature << endl;
403 (*ostr) << "MiuI = " << fMiuI << endl;
404 (*ostr) << "MiuS = " << fMiuS << endl;
405 (*ostr) << "MiuB = " << fMiuB << endl;
406 (*ostr) << "AlphaRange = " << fAlfaRange << endl;
407 (*ostr) << "RapidityRange = " << fRapRange << endl;
408 (*ostr) << "NumberOfIntegrateSamples = 1000000" << endl;
411 void AliGenTherminator::SetModel(const char *model)
413 // Set the freeze-out model to use
414 fFreezeOutModel = model;
415 AliWarning(Form("Selected model %s", fFreezeOutModel.Data()));
416 AliWarning(Form("FOHSLocation is %s", fFOHSlocation.Data()));
419 void AliGenTherminator::SetLhyquidSet(const char *set)
421 // Select one of pregenerated Lhyquid hypersurfaces
422 const char *aroot = gSystem->Getenv("ALICE_ROOT");
423 if (strstr(set, "LHC500C0005")) {
424 AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
425 AliWarning(Form(" Pb-Pb collisions, centrality 0-5 percent"));
426 AliWarning(Form(" initial temperature at tau=1 fm in the center Ti=500 MeV"));
427 AliWarning(Form(" freeze-out criteria Tf=145 MeV"));
428 AliWarning(Form(" for details see $(ALICE_ROOT)/TTherminator/data/LHC500C0005/FO.txt"));
429 fFOHSlocation.Append(aroot);
430 fFOHSlocation.Append("/TTherminator/data/LHC500C0005");
432 if (strstr(set, "LHC500C0510")) {
433 AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
434 AliWarning(Form(" Pb-Pb collisions, centrality 5-10 percent"));
435 AliWarning(Form(" initial temperature at tau=1 fm in the center Ti=500 MeV"));
436 AliWarning(Form(" freeze-out criteria Tf=145 MeV"));
437 AliWarning(Form(" for details see $(ALICE_ROOT)/TTherminator/data/LHC500C0510/FO.txt"));
438 fFOHSlocation.Append(aroot);
439 fFOHSlocation.Append("/TTherminator/data/LHC500C0510");
441 if (strstr(set, "LHC500C1020")) {
442 AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
443 AliWarning(Form(" Pb-Pb collisions, centrality 10-20 percent"));
444 AliWarning(Form(" initial temperature at tau=1 fm in the center Ti=500 MeV"));
445 AliWarning(Form(" freeze-out criteria Tf=145 MeV"));
446 AliWarning(Form(" for details see $(ALICE_ROOT)/TTherminator/data/LHC500C1020/FO.txt"));
447 fFOHSlocation.Append(aroot);
448 fFOHSlocation.Append("/TTherminator/data/LHC500C1020");
450 else if (strstr(set, "LHC500C2030")) {
451 AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
452 AliWarning(Form(" Pb-Pb collisions, centrality 20-30 percent"));
453 AliWarning(Form(" initial temperature at tau=1 fm in the center Ti=500 MeV"));
454 AliWarning(Form(" freeze-out criteria Tf=145 MeV"));
455 AliWarning(Form(" for details see $(ALICE_ROOT)/TTherminator/data/LHC500C2030/FO.txt"));
456 fFOHSlocation.Append(aroot);
457 fFOHSlocation.Append("/TTherminator/data/LHC500C2030");
459 else if (strstr(set, "LHC500C3040")) {
460 AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
461 AliWarning(Form(" Pb-Pb collisions, centrality 30-40 percent"));
462 AliWarning(Form(" initial temperature at tau=1 fm in the center Ti=500 MeV"));
463 AliWarning(Form(" freeze-out criteria Tf=145 MeV"));
464 AliWarning(Form(" for details see $(ALICE_ROOT)/TTherminator/data/LHC500C3040/FO.txt"));
465 fFOHSlocation.Append(aroot);
466 fFOHSlocation.Append("/TTherminator/data/LHC500C3040");
468 else if (strstr(set, "LHC500C4050")) {
469 AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
470 AliWarning(Form(" Pb-Pb collisions, centrality 40-50 percent"));
471 AliWarning(Form(" initial temperature at tau=1 fm in the center Ti=500 MeV"));
472 AliWarning(Form(" freeze-out criteria Tf=145 MeV"));
473 AliWarning(Form(" for details see $(ALICE_ROOT)/TTherminator/data/LHC500C4050/FO.txt"));
474 fFOHSlocation.Append(aroot);
475 fFOHSlocation.Append("/TTherminator/data/LHC500C4050");
477 else if (strstr(set, "LHC276TC0005")) {
478 AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
479 AliWarning(Form(" Pb-Pb collisions, 2.76TeV, centrality 00-05 percent"));
480 AliWarning(Form(" freeze-out criteria Tf=145 MeV"));
481 AliWarning(Form(" for details see $(ALICE_ROOT)/TTherminator/data/LHC276TC0005/FO.txt"));
482 fFOHSlocation.Append(aroot);
483 fFOHSlocation.Append("/TTherminator/data/LHC276TC0005");
485 else if (strstr(set, "LHC276TC0510")) {
486 AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
487 AliWarning(Form(" Pb-Pb collisions, 2.76TeV, centrality 05-10 percent"));
488 AliWarning(Form(" freeze-out criteria Tf=145 MeV"));
489 AliWarning(Form(" for details see $(ALICE_ROOT)/TTherminator/data/LHC276TC0510/FO.txt"));
490 fFOHSlocation.Append(aroot);
491 fFOHSlocation.Append("/TTherminator/data/LHC276TC0510");
493 else if (strstr(set, "LHC276TC1020")) {
494 AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
495 AliWarning(Form(" Pb-Pb collisions, 2.76TeV, centrality 10-20 percent"));
496 AliWarning(Form(" freeze-out criteria Tf=145 MeV"));
497 AliWarning(Form(" for details see $(ALICE_ROOT)/TTherminator/data/LHC276TC1020/FO.txt"));
498 fFOHSlocation.Append(aroot);
499 fFOHSlocation.Append("/TTherminator/data/LHC276TC1020");
501 else if (strstr(set, "LHC276TC2030")) {
502 AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
503 AliWarning(Form(" Pb-Pb collisions, 2.76TeV, centrality 20-30 percent"));
504 AliWarning(Form(" freeze-out criteria Tf=145 MeV"));
505 AliWarning(Form(" for details see $(ALICE_ROOT)/TTherminator/data/LHC276TC2030/FO.txt"));
506 fFOHSlocation.Append(aroot);
507 fFOHSlocation.Append("/TTherminator/data/LHC276TC2030");
509 else if (strstr(set, "LHC276TC3040")) {
510 AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
511 AliWarning(Form(" Pb-Pb collisions, 2.76TeV, centrality 30-40 percent"));
512 AliWarning(Form(" freeze-out criteria Tf=145 MeV"));
513 AliWarning(Form(" for details see $(ALICE_ROOT)/TTherminator/data/LHC276TC3040/FO.txt"));
514 fFOHSlocation.Append(aroot);
515 fFOHSlocation.Append("/TTherminator/data/LHC276TC3040");
517 else if (strstr(set, "LHC276TC4050")) {
518 AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
519 AliWarning(Form(" Pb-Pb collisions, 2.76TeV, centrality 40-50 percent"));
520 AliWarning(Form(" freeze-out criteria Tf=145 MeV"));
521 AliWarning(Form(" for details see $(ALICE_ROOT)/TTherminator/data/LHC276TC4050/FO.txt"));
522 fFOHSlocation.Append(aroot);
523 fFOHSlocation.Append("/TTherminator/data/LHC276TC4050");
526 AliWarning(Form("Did not find Lhyquid set %s", set));
527 AliWarning(Form("Reverting to default: current directory"));
533 void AliGenTherminator::SetLhyquidInputDir(const char *inputdir)
535 // Select Your own Lhyquid hypersurface
536 fFOHSlocation = inputdir;