Update master to aliroot
[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   fFOHSlocation(""),
35   fTemperature(0.1656),
36   fMiuI(-0.0009),
37   fMiuS(0.0),
38   fMiuB(0.0008),
39   fAlfaRange(6.0),
40   fRapRange(4.0),
41   fRhoMax(8.0),
42   fTau(8.0),
43   fBWA(0.0),
44   fBWVt(1.41),
45   fBWDelay(0.0)
46 {
47   // Default constructor
48   fFOHSlocation = "";
49
50   fEnergyCMS = 5500.;
51   fAProjectile = 208;
52   fZProjectile = 82;
53   fProjectile = "A";
54   fATarget = 208;
55   fZTarget = 82;
56   fTarget = "A";
57 }
58 AliGenTherminator::AliGenTherminator(Int_t npart):
59   AliGenMC(npart),
60   fNt(0),
61   fEventNumber(0),
62   fFileName(""),
63   fFreezeOutModel(""),
64   fFOHSlocation(""),
65   fTemperature(0.1656),
66   fMiuI(-0.0009),
67   fMiuS(0.0),
68   fMiuB(0.0008),
69   fAlfaRange(6.0),
70   fRapRange(4.0),
71   fRhoMax(8.0),
72   fTau(8.0),
73   fBWA(0.0),
74   fBWVt(1.41),
75   fBWDelay(0.0)
76 {
77   // Constructor specifying the size of the particle table
78   fNprimaries = 0;
79   fEnergyCMS = 5500.;
80   fAProjectile = 208;
81   fZProjectile = 82;
82   fProjectile = "A";
83   fATarget = 208;
84   fZTarget = 82;
85   fTarget = "A";
86 }
87
88 AliGenTherminator::~AliGenTherminator()
89 {
90   //  AliGenMC::~AliGenMC();
91   //  if (fTherminator) delete fTherminator;
92 }
93
94 void AliGenTherminator::Generate()
95 {
96   // Run single event generation with the Therminator model
97   AliWarning("Generating event from AliGenTherminator");
98
99   Float_t polar[3]    =   {0,0,0};
100   Float_t origin[3]   =   {0,0,0};
101   Float_t time        =   0.;
102   Float_t origin0[3]  =   {0,0,0};
103   Float_t time0       =   0.;
104   Float_t p[3];
105   Float_t mass, energy;
106
107   Int_t nt  = 0;
108   Int_t j, kf, ks, imo;
109   kf = 0;
110
111   Vertex();
112   for (j=0; j < 3; j++) origin0[j] = fVertex[j];
113   time0 = fTime;
114
115   // Generate one event
116
117   ((TTherminator *) fMCEvGen)->GenerateEvent();
118   AliWarning("Generated");
119   ((TTherminator *) fMCEvGen)->ImportParticles(&fParticles);
120
121   Int_t np = fParticles.GetEntriesFast();
122   AliWarning(Form("Imported %d particles", np));
123
124   Int_t *idsOnStack;
125   idsOnStack = new Int_t[np];
126
127   TParticle *iparticle;
128   Double_t evrot = gRandom->Rndm()*TMath::Pi();
129   
130   for (int i = 0; i < np; i++) {
131     iparticle = (TParticle *) fParticles.At(i);
132     Bool_t  hasMother   = (iparticle->GetFirstMother()     >=0);
133     Bool_t  hasDaughter = (iparticle->GetFirstDaughter()   >=0);
134     
135     if (hasDaughter) {
136       // This particle has decayed
137       // It will not be tracked
138       // Add it only once with coorduinates not
139       // smeared with primary vertex position
140
141       kf   = iparticle->GetPdgCode();
142       ks   = iparticle->GetStatusCode();
143       Double_t aphi = TMath::ATan2(iparticle->Py(), iparticle->Px());
144       Double_t arho = TMath::Hypot(iparticle->Px(), iparticle->Py());
145       p[0] = arho*TMath::Cos(aphi + evrot);
146       p[1] = arho*TMath::Sin(aphi + evrot);
147       //     p[0] = iparticle->Px();
148       //     p[1] = iparticle->Py();
149       p[2] = iparticle->Pz();
150       mass = TDatabasePDG::Instance()->GetParticle(kf)->Mass();
151       energy = sqrt(mass*mass + p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
152       
153       Double_t vphi = TMath::ATan2(iparticle->Vy(), iparticle->Vx());
154       Double_t vrho = TMath::Hypot(iparticle->Vx(), iparticle->Vy());
155       origin[0] = vrho*TMath::Cos(vphi + evrot);
156       origin[1] = vrho*TMath::Sin(vphi + evrot);
157       origin[2] = iparticle->Vz();
158       time      = iparticle->T();
159
160       imo = -1;
161       //      TParticle* mother = 0;
162       if (hasMother) {
163         imo = iparticle->GetFirstMother();
164         //      mother = (TParticle *) fParticles.At(imo);
165       } // if has mother   
166       Bool_t tFlag = (!hasDaughter);
167       
168       //      printf("Pushing Track %d with status %d mother %d\n", kf, tFlag, imo>=0?idsOnStack[imo]:imo);
169       PushTrack(tFlag,imo>=0?idsOnStack[imo]:imo,kf,
170                 p[0],p[1],p[2],energy,
171                 origin[0],origin[1],origin[2],time,
172                 polar[0],polar[1],polar[2],
173                 hasMother ? kPDecay:kPNoProcess,nt);
174       idsOnStack[i] = nt;
175       fNprimaries++;
176       KeepTrack(nt);
177     }
178     else {
179       // This is a final state particle
180       // It will be tracked
181       // Add it TWICE to the stack !!!
182       // First time with event-wide coordicates (for femtoscopy) - 
183       //   this one will not be tracked
184       // Second time with event-wide ccordiantes and vertex smearing
185       //   this one will be tracked
186       
187       kf   = iparticle->GetPdgCode();
188       ks   = iparticle->GetStatusCode();
189       Double_t aphi = TMath::ATan2(iparticle->Py(), iparticle->Px());
190       Double_t arho = TMath::Hypot(iparticle->Px(), iparticle->Py());
191       p[0] = arho*TMath::Cos(aphi + evrot);
192       p[1] = arho*TMath::Sin(aphi + evrot);
193       //     p[0] = iparticle->Px();
194       //     p[1] = iparticle->Py();
195       p[2] = iparticle->Pz();
196       mass = TDatabasePDG::Instance()->GetParticle(kf)->Mass();
197       energy = sqrt(mass*mass + p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
198       
199       Double_t vphi = TMath::ATan2(iparticle->Vy(), iparticle->Vx());
200       Double_t vrho = TMath::Hypot(iparticle->Vx(), iparticle->Vy());
201       origin[0] = vrho*TMath::Cos(vphi + evrot);
202       origin[1] = vrho*TMath::Sin(vphi + evrot);
203       origin[2] = iparticle->Vz();
204       time      = iparticle->T();
205
206       imo = -1;
207       //      TParticle* mother = 0;
208       if (hasMother) {
209         imo = iparticle->GetFirstMother();
210         //      mother = (TParticle *) fParticles.At(imo);
211       } // if has mother   
212       Bool_t tFlag = (hasDaughter);
213       
214 //       printf("Found mother %i with true id %i\n", imo, imo>=0?idsOnStack[imo]:imo);
215 //       printf("Pushing Track %d with status %d mother %d\n", kf, tFlag, imo>=0?idsOnStack[imo]:imo);
216       PushTrack(tFlag,imo>=0?idsOnStack[imo]:imo,kf,
217                 p[0],p[1],p[2],energy,
218                 origin[0],origin[1],origin[2],time,
219                 polar[0],polar[1],polar[2],
220                 hasMother ? kPDecay:kPNoProcess,nt);
221       idsOnStack[i] = nt;
222       fNprimaries++;
223       KeepTrack(nt);
224
225       origin[0] = origin0[0]+vrho*TMath::Cos(vphi + evrot);
226       origin[1] = origin0[1]+vrho*TMath::Sin(vphi + evrot);
227       origin[2] = origin0[2]+iparticle->Vz();
228       time      = time0+iparticle->T();
229
230       imo = nt;
231       //      mother = (TParticle *) fParticles.At(nt);
232       tFlag = (!hasDaughter);
233       
234 //       printf("Pushing Track %d with status %d mother %d\n", kf, tFlag, imo);
235       PushTrack(tFlag,imo,kf,
236                 p[0],p[1],p[2],energy,
237                 origin[0],origin[1],origin[2],time,
238                 polar[0],polar[1],polar[2],
239                 hasMother ? kPDecay:kPNoProcess,nt);
240       fNprimaries++;
241       KeepTrack(nt);
242     }
243   }
244
245
246
247   SetHighWaterMark(fNprimaries);
248
249   TArrayF eventVertex;
250   eventVertex.Set(3);
251   eventVertex[0] = origin0[0];
252   eventVertex[1] = origin0[1];
253   eventVertex[2] = origin0[2];
254   Float_t eventTime = time0;
255
256 // Builds the event header, to be called after each event
257   AliGenEventHeader* header = new AliGenHijingEventHeader("Therminator");
258
259   // Header
260   //  AliGenEventHeader* header = new AliGenEventHeader("Therminator");
261   // Event Vertex
262 //   header->SetPrimaryVertex(eventVertex);
263 //   header->SetNProduced(fNprimaries);
264
265   ((AliGenHijingEventHeader*) header)->SetNProduced(fNprimaries);
266   ((AliGenHijingEventHeader*) header)->SetPrimaryVertex(eventVertex);
267   ((AliGenHijingEventHeader*) header)->SetInteractionTime(eventTime);
268   ((AliGenHijingEventHeader*) header)->SetImpactParameter(0.0);
269   ((AliGenHijingEventHeader*) header)->SetTotalEnergy(0.0);
270   ((AliGenHijingEventHeader*) header)->SetHardScatters(0);
271   ((AliGenHijingEventHeader*) header)->SetParticipants(0, 0);
272   ((AliGenHijingEventHeader*) header)->SetCollisions(0, 0, 0, 0);
273   ((AliGenHijingEventHeader*) header)->SetSpectators(0, 0, 0, 0);
274   ((AliGenHijingEventHeader*) header)->SetReactionPlaneAngle(evrot);
275
276
277 // 4-momentum vectors of the triggered jets.
278 //
279 // Before final state gluon radiation.
280 //     TLorentzVector* jet1 = new TLorentzVector(fHijing->GetHINT1(21), 
281 //                                            fHijing->GetHINT1(22),
282 //                                            fHijing->GetHINT1(23),
283 //                                            fHijing->GetHINT1(24));
284
285 //     TLorentzVector* jet2 = new TLorentzVector(fHijing->GetHINT1(31), 
286 //                                            fHijing->GetHINT1(32),
287 //                                            fHijing->GetHINT1(33),
288 //                                            fHijing->GetHINT1(34));
289 // // After final state gluon radiation.
290 //     TLorentzVector* jet3 = new TLorentzVector(fHijing->GetHINT1(26), 
291 //                                            fHijing->GetHINT1(27),
292 //                                            fHijing->GetHINT1(28),
293 //                                            fHijing->GetHINT1(29));
294
295 //     TLorentzVector* jet4 = new TLorentzVector(fHijing->GetHINT1(36), 
296 //                                            fHijing->GetHINT1(37),
297 //                                            fHijing->GetHINT1(38),
298 //                                            fHijing->GetHINT1(39));
299 //     ((AliGenHijingEventHeader*) header)->SetJets(jet1, jet2, jet3, jet4);
300 // Bookkeeping for kinematic bias
301 //     ((AliGenHijingEventHeader*) header)->SetTrials(fTrials);
302 // Event Vertex
303   header->SetPrimaryVertex(fVertex);
304   header->SetInteractionTime(fTime);
305   AddHeader(header);
306   fCollisionGeometry = (AliGenHijingEventHeader*)  header;
307
308   delete [] idsOnStack;
309
310   //  gAlice->SetGenEventHeader(header); 
311 }
312
313 void AliGenTherminator::Init()
314 {
315   // Initialize global variables and
316   // particle and decay tables
317   if (fFileName.Length() == 0)
318     fFileName = "event.out";
319   ReadShareParticleTable();
320
321   SetMC(new TTherminator());
322   
323   AliGenMC::Init();
324   ((TTherminator *) fMCEvGen)->Initialize();
325 }
326 void AliGenTherminator::SetFileName(const char *infilename)
327 {
328   // Set parameter filename
329   fFileName = infilename;
330 }
331 void AliGenTherminator::SetEventNumberInFile(int evnum)
332 {
333   // Set number of events to generate - default: 1
334   fEventNumber = evnum;
335 }
336
337 void AliGenTherminator::ReadShareParticleTable()
338 {
339   // Read in particle table from share
340   // and add missing particle type to TDatabasePDG
341
342   char str[50];
343   char str1[200];
344     
345   TDatabasePDG *tInstance = TDatabasePDG::Instance();
346   TParticlePDG *tParticleType;
347
348   AliWarning(Form("Reading particle types from particles.data"));
349
350   TString aroot = gSystem->Getenv("ALICE_ROOT");
351   ifstream in((aroot+"/TTherminator/data/SHARE/particles.data").Data());
352   //  ifstream in("particles.data");
353   
354   int charge;
355     
356   int number=0;
357   if ((in) && (in.is_open()))
358     {
359       //START OF HEAD-LINE
360       in.ignore(200,'\n');
361       in.ignore(200,'\n');
362       in.ignore(200,'\n');
363       //END OF HEAD-LINE
364       
365       while (in>>str)
366         {
367           if (/*(*str == '#')||*/(*str<65)||(*str>122))
368             {
369               in.getline(str1,200);
370               continue;
371             }
372           double mass, gamma, spin, tI3, tI, q, s, aq, as, c, ac, mc;
373           
374           in>>mass>>gamma>>spin>>tI>>tI3>>q>>s>>aq>>as>>c>>ac>>mc;
375           number++;
376           tParticleType = tInstance->GetParticle((int) mc);
377           if (!tParticleType) {
378             charge = 0;
379             if (strstr(str, "plu")) charge = 1;
380             if (strstr(str, "min")) charge = -1;
381             if (strstr(str, "plb")) charge = -1;
382             if (strstr(str, "mnb")) charge = 1;
383             if (strstr(str, "plp")) charge = 2;
384             if (strstr(str, "ppb")) charge = -2;
385             tInstance->AddParticle(str, str, mass, gamma == 0.0 ? 1:0, gamma, charge , "meson", (int) mc);
386             AliWarning(Form("Added particle %s with PDG PID %d charge %d", str, (int) mc, charge));
387             //      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));
388
389           }
390         }
391       in.close();
392     }
393   CreateTherminatorInputFile();
394 }
395
396 void AliGenTherminator::CreateTherminatorInputFile()
397 {
398   // Create Therminator input file
399   const char *aroot = gSystem->Getenv("ALICE_ROOT");
400   ofstream *ostr = new ofstream("therminator.in");
401   (*ostr) << "NumberOfEvents = 1" << endl;
402   (*ostr) << "Randomize = 1" << endl;
403   (*ostr) << "TableType = SHARE" << endl;
404   (*ostr) << "InputDirSHARE = "<< aroot << "/TTherminator/data/SHARE" << endl;
405   (*ostr) << "EventOutputFile = " << fFileName.Data() << endl;
406   (*ostr) << "FOHSLocation = " << fFOHSlocation.Data() << endl;
407   (*ostr) << "FreezeOutModel = " << fFreezeOutModel.Data() << endl;
408   (*ostr) << "BWVt = " << fBWVt << endl;
409   (*ostr) << "Tau = " << fTau << endl;
410   (*ostr) << "RhoMax = " << fRhoMax << endl;
411   (*ostr) << "Temperature = " << fTemperature << endl;
412   (*ostr) << "MiuI = " << fMiuI << endl;
413   (*ostr) << "MiuS = " << fMiuS << endl;
414   (*ostr) << "MiuB = " << fMiuB << endl;
415   (*ostr) << "AlphaRange = " << fAlfaRange << endl;
416   (*ostr) << "RapidityRange = " << fRapRange << endl;
417   (*ostr) << "NumberOfIntegrateSamples = 1000000" << endl;
418 }
419
420 void AliGenTherminator::SetModel(const char *model)
421 {
422   // Set the freeze-out model to use
423   fFreezeOutModel = model;
424   AliWarning(Form("Selected model %s", fFreezeOutModel.Data()));
425   AliWarning(Form("FOHSLocation is %s", fFOHSlocation.Data()));
426 }
427
428 void AliGenTherminator::SetLhyquidSet(const char *set)
429 {
430   // Select one of pregenerated Lhyquid hypersurfaces
431   const char *aroot = gSystem->Getenv("ALICE_ROOT");
432   if (strstr(set, "LHC500C0005")) {
433     AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
434     AliWarning(Form("  Pb-Pb collisions, centrality 0-5 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/LHC500C0005/FO.txt"));
438     fFOHSlocation.Append(aroot);
439     fFOHSlocation.Append("/TTherminator/data/LHC500C0005");
440   }
441   if (strstr(set, "LHC500C0510")) {
442     AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
443     AliWarning(Form("  Pb-Pb collisions, centrality 5-10 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/LHC500C0510/FO.txt"));
447     fFOHSlocation.Append(aroot);
448     fFOHSlocation.Append("/TTherminator/data/LHC500C0510");
449   }
450   if (strstr(set, "LHC500C1020")) {
451     AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
452     AliWarning(Form("  Pb-Pb collisions, centrality 10-20 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/LHC500C1020/FO.txt"));
456     fFOHSlocation.Append(aroot);
457     fFOHSlocation.Append("/TTherminator/data/LHC500C1020");
458   }
459   else if (strstr(set, "LHC500C2030")) {
460     AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
461     AliWarning(Form("  Pb-Pb collisions, centrality 20-30 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/LHC500C2030/FO.txt"));
465     fFOHSlocation.Append(aroot);
466     fFOHSlocation.Append("/TTherminator/data/LHC500C2030");
467   }
468   else if (strstr(set, "LHC500C3040")) {
469     AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
470     AliWarning(Form("  Pb-Pb collisions, centrality 30-40 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/LHC500C3040/FO.txt"));
474     fFOHSlocation.Append(aroot);
475     fFOHSlocation.Append("/TTherminator/data/LHC500C3040");
476   }
477   else if (strstr(set, "LHC500C4050")) {
478     AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
479     AliWarning(Form("  Pb-Pb collisions, centrality 40-50 percent"));
480     AliWarning(Form("  initial temperature at tau=1 fm in the center Ti=500 MeV"));
481     AliWarning(Form("  freeze-out criteria Tf=145 MeV"));
482     AliWarning(Form("  for details see $(ALICE_ROOT)/TTherminator/data/LHC500C4050/FO.txt"));
483     fFOHSlocation.Append(aroot);
484     fFOHSlocation.Append("/TTherminator/data/LHC500C4050");
485   }
486   else if (strstr(set, "LHC276TC0005")) {
487     AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
488     AliWarning(Form("  Pb-Pb collisions, 2.76TeV, centrality 00-05 percent"));
489     AliWarning(Form("  freeze-out criteria Tf=145 MeV"));
490     AliWarning(Form("  for details see $(ALICE_ROOT)/TTherminator/data/LHC276TC0005/FO.txt"));
491     fFOHSlocation.Append(aroot);
492     fFOHSlocation.Append("/TTherminator/data/LHC276TC0005");
493   }
494   else if (strstr(set, "LHC276TC0510")) {
495     AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
496     AliWarning(Form("  Pb-Pb collisions, 2.76TeV, centrality 05-10 percent"));
497     AliWarning(Form("  freeze-out criteria Tf=145 MeV"));
498     AliWarning(Form("  for details see $(ALICE_ROOT)/TTherminator/data/LHC276TC0510/FO.txt"));
499     fFOHSlocation.Append(aroot);
500     fFOHSlocation.Append("/TTherminator/data/LHC276TC0510");
501   }
502   else if (strstr(set, "LHC276TC1020")) {
503     AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
504     AliWarning(Form("  Pb-Pb collisions, 2.76TeV, centrality 10-20 percent"));
505     AliWarning(Form("  freeze-out criteria Tf=145 MeV"));
506     AliWarning(Form("  for details see $(ALICE_ROOT)/TTherminator/data/LHC276TC1020/FO.txt"));
507     fFOHSlocation.Append(aroot);
508     fFOHSlocation.Append("/TTherminator/data/LHC276TC1020");
509   }
510   else if (strstr(set, "LHC276TC2030")) {
511     AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
512     AliWarning(Form("  Pb-Pb collisions, 2.76TeV, centrality 20-30 percent"));
513     AliWarning(Form("  freeze-out criteria Tf=145 MeV"));
514     AliWarning(Form("  for details see $(ALICE_ROOT)/TTherminator/data/LHC276TC2030/FO.txt"));
515     fFOHSlocation.Append(aroot);
516     fFOHSlocation.Append("/TTherminator/data/LHC276TC2030");
517   }
518   else if (strstr(set, "LHC276TC3040")) {
519     AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
520     AliWarning(Form("  Pb-Pb collisions, 2.76TeV, centrality 30-40 percent"));
521     AliWarning(Form("  freeze-out criteria Tf=145 MeV"));
522     AliWarning(Form("  for details see $(ALICE_ROOT)/TTherminator/data/LHC276TC3040/FO.txt"));
523     fFOHSlocation.Append(aroot);
524     fFOHSlocation.Append("/TTherminator/data/LHC276TC3040");
525   }
526   else if (strstr(set, "LHC276TC4050")) {
527     AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
528     AliWarning(Form("  Pb-Pb collisions, 2.76TeV, centrality 40-50 percent"));
529     AliWarning(Form("  freeze-out criteria Tf=145 MeV"));
530     AliWarning(Form("  for details see $(ALICE_ROOT)/TTherminator/data/LHC276TC4050/FO.txt"));
531     fFOHSlocation.Append(aroot);
532     fFOHSlocation.Append("/TTherminator/data/LHC276TC4050");
533   }
534   else {
535     AliWarning(Form("Did not find Lhyquid set %s", set));
536     AliWarning(Form("Reverting to default: current directory"));
537     fFOHSlocation += "";
538   }
539
540 }
541
542 void AliGenTherminator::SetLhyquidInputDir(const char *inputdir)
543 {
544   // Select Your own Lhyquid hypersurface
545   fFOHSlocation = inputdir;
546 }