Completing the setup for TTherminator model.
[u/mrichter/AliRoot.git] / TTherminator / AliGenTherminator.cxx
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
4 // It has an option to use the Lhyquid3D input freeze-out
5 // hypersurface
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"
19 #include "AliGenHijingEventHeader.h"
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;
106   Double_t evrot = gRandom->Rndm()*TMath::Pi();
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();
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();
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
155 // Builds the event header, to be called after each event
156   AliGenEventHeader* header = new AliGenHijingEventHeader("Therminator");
157
158   // Header
159   //  AliGenEventHeader* header = new AliGenEventHeader("Therminator");
160   // Event Vertex
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
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"));
244
245   TString aroot = gSystem->Getenv("ALICE_ROOT");
246   ifstream in((aroot+"/TTherminator/data/SHARE/particles.data").Data());
247   //  ifstream in("particles.data");
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
294   const char *aroot = gSystem->Getenv("ALICE_ROOT");
295   ofstream *ostr = new ofstream("therminator.in");
296   (*ostr) << "NumberOfEvents = 1" << endl;
297   (*ostr) << "Randomize = 1" << endl;
298   (*ostr) << "TableType = SHARE" << endl;
299   (*ostr) << "InputDirSHARE = "<< aroot << "/TTherminator/data/SHARE" << endl;
300   (*ostr) << "EventOutputFile = " << fFileName.Data() << endl;
301   (*ostr) << "FOHSLocation = " << fFOHSlocation.Data() << endl;
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 }
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 }