Completing the setup for TTherminator model.
[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
47 fParticles = new TClonesArray("TParticle",20000);
48}
49AliGenTherminator::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
72AliGenTherminator::~AliGenTherminator()
73{
74 // AliGenMC::~AliGenMC();
75 // if (fTherminator) delete fTherminator;
76}
77
78void 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
208void 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}
221void AliGenTherminator::SetFileName(const char *infilename)
222{
223 // Set parameter filename
224 fFileName = infilename;
225}
226void AliGenTherminator::SetEventNumberInFile(int evnum)
227{
228 // Set number of events to generate - default: 1
229 fEventNumber = evnum;
230}
231
232void 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
291void 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
315void AliGenTherminator::SetModel(const char *model)
316{
317 // Set the freeze-out model to use
318 fFreezeOutModel = model;
319}
ae89e57a 320
321void 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
350void AliGenTherminator::SetLhyquidInputDir(const char *inputdir)
351{
352 // Select Your own Lhyquid hypersurface
353 fFOHSlocation = inputdir;
354}