Adding Optional Output
[u/mrichter/AliRoot.git] / TTherminator / AliGenTherminator.cxx
CommitLineData
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
24ClassImp(AliGenTherminator)
25
26using namespace std;
27
28AliGenTherminator::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}
48AliGenTherminator::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
70AliGenTherminator::~AliGenTherminator()
71{
72 // AliGenMC::~AliGenMC();
73 // if (fTherminator) delete fTherminator;
74}
75
76void 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
146 printf("Pushing Track %d with status %d mother %d\n", kf, tFlag, imo);
147 PushTrack(tFlag,imo,kf,
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
deb008f3 191 printf("Found mother %i with true id %i\n", imo, idsOnStack[imo]);
192 printf("Pushing Track %d with status %d mother %d\n", kf, tFlag, idsOnStack[imo]);
c9a371fc 193 PushTrack(tFlag,imo,kf,
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
286void 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}
299void AliGenTherminator::SetFileName(const char *infilename)
300{
301 // Set parameter filename
302 fFileName = infilename;
303}
304void AliGenTherminator::SetEventNumberInFile(int evnum)
305{
306 // Set number of events to generate - default: 1
307 fEventNumber = evnum;
308}
309
310void 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
369void 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
393void AliGenTherminator::SetModel(const char *model)
394{
395 // Set the freeze-out model to use
396 fFreezeOutModel = model;
397}
ae89e57a 398
399void 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"));
405 AliWarning(Form(" Pb-Pb collisions, centrality 0-5%"));
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"));
414 AliWarning(Form(" Pb-Pb collisions, centrality 20-30%"));
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
428void AliGenTherminator::SetLhyquidInputDir(const char *inputdir)
429{
430 // Select Your own Lhyquid hypersurface
431 fFOHSlocation = inputdir;
432}