[u/mrichter/AliRoot.git] / TTherminator / AliGenTherminator.cxx

// ALICE event generator based on the THERMINATOR model
// It reads the test output of the model and puts it onto
// the stack
// It has an option to use the Lhyquid3D input freeze-out
// hypersurface
// Author: Adam.Kisiel@cern.ch

#include <iostream>
#include <fstream>
#include <sstream>
#include <TClonesArray.h>
#include <TMCProcess.h>
#include <TDatabasePDG.h>
#include <TParticle.h>

#include "AliConst.h"
#include "AliDecayer.h"
#include "AliGenEventHeader.h"
#include "AliGenHijingEventHeader.h"
#include "AliGenTherminator.h"
#include "AliLog.h"
#include "AliRun.h"

ClassImp(AliGenTherminator)

using namespace std;

AliGenTherminator::AliGenTherminator():
  AliGenMC(),
  fNt(0),
  fEventNumber(0),
  fFileName(""),
  fFreezeOutModel(""),
  fFOHSlocation(""),
  fTemperature(0.1656),
  fMiuI(-0.0009),
  fMiuS(0.0),
  fMiuB(0.0008),
  fAlfaRange(6.0),
  fRapRange(4.0),
  fRhoMax(8.0),
  fTau(8.0),
  fBWA(0.0),
  fBWVt(1.41),
  fBWDelay(0.0)
{
  // Default constructor
  fFOHSlocation = "";

  fEnergyCMS = 5500.;
  fAProjectile = 208;
  fZProjectile = 82;
  fProjectile = "A";
  fATarget = 208;
  fZTarget = 82;
  fTarget = "A";
}
AliGenTherminator::AliGenTherminator(Int_t npart):
  AliGenMC(npart),
  fNt(0),
  fEventNumber(0),
  fFileName(""),
  fFreezeOutModel(""),
  fFOHSlocation(""),
  fTemperature(0.1656),
  fMiuI(-0.0009),
  fMiuS(0.0),
  fMiuB(0.0008),
  fAlfaRange(6.0),
  fRapRange(4.0),
  fRhoMax(8.0),
  fTau(8.0),
  fBWA(0.0),
  fBWVt(1.41),
  fBWDelay(0.0)
{
  // Constructor specifying the size of the particle table
  fNprimaries = 0;
  fEnergyCMS = 5500.;
  fAProjectile = 208;
  fZProjectile = 82;
  fProjectile = "A";
  fATarget = 208;
  fZTarget = 82;
  fTarget = "A";
}

AliGenTherminator::~AliGenTherminator()
{
  //  AliGenMC::~AliGenMC();
  //  if (fTherminator) delete fTherminator;
}

void AliGenTherminator::Generate()
{
  // Run single event generation with the Therminator model
  AliWarning("Generating event from AliGenTherminator");

  Float_t polar[3]    =   {0,0,0};
  Float_t origin[3]   =   {0,0,0};
  Float_t time        =   0.;
  Float_t origin0[3]  =   {0,0,0};
  Float_t time0       =   0.;
  Float_t p[3];
  Float_t mass, energy;

  Int_t nt  = 0;
  Int_t j, kf, ks, imo;
  kf = 0;

  Vertex();
  for (j=0; j < 3; j++) origin0[j] = fVertex[j];
  time0 = fTime;

  // Generate one event

  ((TTherminator *) fMCEvGen)->GenerateEvent();
  AliWarning("Generated");
  ((TTherminator *) fMCEvGen)->ImportParticles(&fParticles);

  Int_t np = fParticles.GetEntriesFast();
  AliWarning(Form("Imported %d particles", np));

  Int_t *idsOnStack;
  idsOnStack = new Int_t[np];

  TParticle *iparticle;
  Double_t evrot = gRandom->Rndm()*TMath::Pi();
  
  for (int i = 0; i < np; i++) {
    iparticle = (TParticle *) fParticles.At(i);
    Bool_t  hasMother   = (iparticle->GetFirstMother()     >=0);
    Bool_t  hasDaughter = (iparticle->GetFirstDaughter()   >=0);
    
    if (hasDaughter) {
      // This particle has decayed
      // It will not be tracked
      // Add it only once with coorduinates not
      // smeared with primary vertex position

      kf   = iparticle->GetPdgCode();
      ks   = iparticle->GetStatusCode();
      Double_t aphi = TMath::ATan2(iparticle->Py(), iparticle->Px());
      Double_t arho = TMath::Hypot(iparticle->Px(), iparticle->Py());
      p[0] = arho*TMath::Cos(aphi + evrot);
      p[1] = arho*TMath::Sin(aphi + evrot);
      //     p[0] = iparticle->Px();
      //     p[1] = iparticle->Py();
      p[2] = iparticle->Pz();
      mass = TDatabasePDG::Instance()->GetParticle(kf)->Mass();
      energy = sqrt(mass*mass + p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
      
      Double_t vphi = TMath::ATan2(iparticle->Vy(), iparticle->Vx());
      Double_t vrho = TMath::Hypot(iparticle->Vx(), iparticle->Vy());
      origin[0] = vrho*TMath::Cos(vphi + evrot);
      origin[1] = vrho*TMath::Sin(vphi + evrot);
      origin[2] = iparticle->Vz();
      time      = iparticle->T();

      imo = -1;
      //      TParticle* mother = 0;
      if (hasMother) {
	imo = iparticle->GetFirstMother();
	//	mother = (TParticle *) fParticles.At(imo);
      } // if has mother   
      Bool_t tFlag = (!hasDaughter);
      
      //      printf("Pushing Track %d with status %d mother %d\n", kf, tFlag, imo>=0?idsOnStack[imo]:imo);
      PushTrack(tFlag,imo>=0?idsOnStack[imo]:imo,kf,
		p[0],p[1],p[2],energy,
		origin[0],origin[1],origin[2],time,
		polar[0],polar[1],polar[2],
		hasMother ? kPDecay:kPNoProcess,nt);
      idsOnStack[i] = nt;
      fNprimaries++;
      KeepTrack(nt);
    }
    else {
      // This is a final state particle
      // It will be tracked
      // Add it TWICE to the stack !!!
      // First time with event-wide coordicates (for femtoscopy) - 
      //   this one will not be tracked
      // Second time with event-wide ccordiantes and vertex smearing
      //   this one will be tracked
      
      kf   = iparticle->GetPdgCode();
      ks   = iparticle->GetStatusCode();
      Double_t aphi = TMath::ATan2(iparticle->Py(), iparticle->Px());
      Double_t arho = TMath::Hypot(iparticle->Px(), iparticle->Py());
      p[0] = arho*TMath::Cos(aphi + evrot);
      p[1] = arho*TMath::Sin(aphi + evrot);
      //     p[0] = iparticle->Px();
      //     p[1] = iparticle->Py();
      p[2] = iparticle->Pz();
      mass = TDatabasePDG::Instance()->GetParticle(kf)->Mass();
      energy = sqrt(mass*mass + p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
      
      Double_t vphi = TMath::ATan2(iparticle->Vy(), iparticle->Vx());
      Double_t vrho = TMath::Hypot(iparticle->Vx(), iparticle->Vy());
      origin[0] = vrho*TMath::Cos(vphi + evrot);
      origin[1] = vrho*TMath::Sin(vphi + evrot);
      origin[2] = iparticle->Vz();
      time      = iparticle->T();

      imo = -1;
      //      TParticle* mother = 0;
      if (hasMother) {
	imo = iparticle->GetFirstMother();
	//	mother = (TParticle *) fParticles.At(imo);
      } // if has mother   
      Bool_t tFlag = (hasDaughter);
      
//       printf("Found mother %i with true id %i\n", imo, imo>=0?idsOnStack[imo]:imo);
//       printf("Pushing Track %d with status %d mother %d\n", kf, tFlag, imo>=0?idsOnStack[imo]:imo);
      PushTrack(tFlag,imo>=0?idsOnStack[imo]:imo,kf,
		p[0],p[1],p[2],energy,
		origin[0],origin[1],origin[2],time,
		polar[0],polar[1],polar[2],
		hasMother ? kPDecay:kPNoProcess,nt);
      idsOnStack[i] = nt;
      fNprimaries++;
      KeepTrack(nt);

      origin[0] = origin0[0]+vrho*TMath::Cos(vphi + evrot);
      origin[1] = origin0[1]+vrho*TMath::Sin(vphi + evrot);
      origin[2] = origin0[2]+iparticle->Vz();
      time      = time0+iparticle->T();

      imo = nt;
      //      mother = (TParticle *) fParticles.At(nt);
      tFlag = (!hasDaughter);
      
//       printf("Pushing Track %d with status %d mother %d\n", kf, tFlag, imo);
      PushTrack(tFlag,imo,kf,
		p[0],p[1],p[2],energy,
		origin[0],origin[1],origin[2],time,
		polar[0],polar[1],polar[2],
		hasMother ? kPDecay:kPNoProcess,nt);
      fNprimaries++;
      KeepTrack(nt);
    }
  }


  SetHighWaterMark(fNprimaries);

  TArrayF eventVertex;
  eventVertex.Set(3);
  eventVertex[0] = origin0[0];
  eventVertex[1] = origin0[1];
  eventVertex[2] = origin0[2];
  Float_t eventTime = time0;

// Builds the event header, to be called after each event
  AliGenEventHeader* header = new AliGenHijingEventHeader("Therminator");

  // Header
  //  AliGenEventHeader* header = new AliGenEventHeader("Therminator");
  // Event Vertex
//   header->SetPrimaryVertex(eventVertex);
//   header->SetNProduced(fNprimaries);

  ((AliGenHijingEventHeader*) header)->SetNProduced(fNprimaries);
  ((AliGenHijingEventHeader*) header)->SetPrimaryVertex(eventVertex);
  ((AliGenHijingEventHeader*) header)->SetInteractionTime(eventTime);
  ((AliGenHijingEventHeader*) header)->SetImpactParameter(0.0);
  ((AliGenHijingEventHeader*) header)->SetTotalEnergy(0.0);
  ((AliGenHijingEventHeader*) header)->SetHardScatters(0);
  ((AliGenHijingEventHeader*) header)->SetParticipants(0, 0);
  ((AliGenHijingEventHeader*) header)->SetCollisions(0, 0, 0, 0);
  ((AliGenHijingEventHeader*) header)->SetSpectators(0, 0, 0, 0);
  ((AliGenHijingEventHeader*) header)->SetReactionPlaneAngle(evrot);


// 4-momentum vectors of the triggered jets.
//
// Before final state gluon radiation.
//     TLorentzVector* jet1 = new TLorentzVector(fHijing->GetHINT1(21), 
// 					      fHijing->GetHINT1(22),
// 					      fHijing->GetHINT1(23),
// 					      fHijing->GetHINT1(24));

//     TLorentzVector* jet2 = new TLorentzVector(fHijing->GetHINT1(31), 
// 					      fHijing->GetHINT1(32),
// 					      fHijing->GetHINT1(33),
// 					      fHijing->GetHINT1(34));
// // After final state gluon radiation.
//     TLorentzVector* jet3 = new TLorentzVector(fHijing->GetHINT1(26), 
// 					      fHijing->GetHINT1(27),
// 					      fHijing->GetHINT1(28),
// 					      fHijing->GetHINT1(29));

//     TLorentzVector* jet4 = new TLorentzVector(fHijing->GetHINT1(36), 
// 					      fHijing->GetHINT1(37),
// 					      fHijing->GetHINT1(38),
// 					      fHijing->GetHINT1(39));
//     ((AliGenHijingEventHeader*) header)->SetJets(jet1, jet2, jet3, jet4);
// Bookkeeping for kinematic bias
//     ((AliGenHijingEventHeader*) header)->SetTrials(fTrials);
// Event Vertex
  header->SetPrimaryVertex(fVertex);
  header->SetInteractionTime(fTime);
  AddHeader(header);
  fCollisionGeometry = (AliGenHijingEventHeader*)  header;

  delete [] idsOnStack;

  //  gAlice->SetGenEventHeader(header); 
}

void AliGenTherminator::Init()
{
  // Initialize global variables and
  // particle and decay tables
  if (fFileName.Length() == 0)
    fFileName = "event.out";
  ReadShareParticleTable();

  SetMC(new TTherminator());
  
  AliGenMC::Init();
  ((TTherminator *) fMCEvGen)->Initialize();
}
void AliGenTherminator::SetFileName(const char *infilename)
{
  // Set parameter filename
  fFileName = infilename;
}
void AliGenTherminator::SetEventNumberInFile(int evnum)
{
  // Set number of events to generate - default: 1
  fEventNumber = evnum;
}

void AliGenTherminator::ReadShareParticleTable()
{
  // Read in particle table from share
  // and add missing particle type to TDatabasePDG

  char str[50];
  char str1[200];
    
  TDatabasePDG *tInstance = TDatabasePDG::Instance();
  TParticlePDG *tParticleType;

  AliWarning(Form("Reading particle types from particles.data"));

  TString aroot = gSystem->Getenv("ALICE_ROOT");
  ifstream in((aroot+"/TTherminator/data/SHARE/particles.data").Data());
  //  ifstream in("particles.data");
  
  int charge;
    
  int number=0;
  if ((in) && (in.is_open()))
    {
      //START OF HEAD-LINE
      in.ignore(200,'\n');
      in.ignore(200,'\n');
      in.ignore(200,'\n');
      //END OF HEAD-LINE
      
      while (in>>str)
	{
	  if (/*(*str == '#')||*/(*str<65)||(*str>122))
	    {
	      in.getline(str1,200);
	      continue;
	    }
	  double mass, gamma, spin, tI3, tI, q, s, aq, as, c, ac, mc;
	  
	  in>>mass>>gamma>>spin>>tI>>tI3>>q>>s>>aq>>as>>c>>ac>>mc;
	  number++;
	  tParticleType = tInstance->GetParticle((int) mc);
	  if (!tParticleType) {
	    charge = 0;
	    if (strstr(str, "plu")) charge = 1;
	    if (strstr(str, "min")) charge = -1;
	    if (strstr(str, "plb")) charge = -1;
	    if (strstr(str, "mnb")) charge = 1;
	    if (strstr(str, "plp")) charge = 2;
	    if (strstr(str, "ppb")) charge = -2;
	    tInstance->AddParticle(str, str, mass, gamma == 0.0 ? 1:0, gamma, charge , "meson", (int) mc);
	    AliWarning(Form("Added particle %s with PDG PID %d charge %d", str, (int) mc, charge));
	    //	    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));

	  }
	}
      in.close();
    }
  CreateTherminatorInputFile();
}

void AliGenTherminator::CreateTherminatorInputFile()
{
  // Create Therminator input file
  const char *aroot = gSystem->Getenv("ALICE_ROOT");
  ofstream *ostr = new ofstream("therminator.in");
  (*ostr) << "NumberOfEvents = 1" << endl;
  (*ostr) << "Randomize = 1" << endl;
  (*ostr) << "TableType = SHARE" << endl;
  (*ostr) << "InputDirSHARE = "<< aroot << "/TTherminator/data/SHARE" << endl;
  (*ostr) << "EventOutputFile = " << fFileName.Data() << endl;
  (*ostr) << "FOHSLocation = " << fFOHSlocation.Data() << endl;
  (*ostr) << "FreezeOutModel = " << fFreezeOutModel.Data() << endl;
  (*ostr) << "BWVt = " << fBWVt << endl;
  (*ostr) << "Tau = " << fTau << endl;
  (*ostr) << "RhoMax = " << fRhoMax << endl;
  (*ostr) << "Temperature = " << fTemperature << endl;
  (*ostr) << "MiuI = " << fMiuI << endl;
  (*ostr) << "MiuS = " << fMiuS << endl;
  (*ostr) << "MiuB = " << fMiuB << endl;
  (*ostr) << "AlphaRange = " << fAlfaRange << endl;
  (*ostr) << "RapidityRange = " << fRapRange << endl;
  (*ostr) << "NumberOfIntegrateSamples = 1000000" << endl;
}

void AliGenTherminator::SetModel(const char *model)
{
  // Set the freeze-out model to use
  fFreezeOutModel = model;
  AliWarning(Form("Selected model %s", fFreezeOutModel.Data()));
  AliWarning(Form("FOHSLocation is %s", fFOHSlocation.Data()));
}

void AliGenTherminator::SetLhyquidSet(const char *set)
{
  // Select one of pregenerated Lhyquid hypersurfaces
  const char *aroot = gSystem->Getenv("ALICE_ROOT");
  if (strstr(set, "LHC500C0005")) {
    AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
    AliWarning(Form("  Pb-Pb collisions, centrality 0-5 percent"));
    AliWarning(Form("  initial temperature at tau=1 fm in the center Ti=500 MeV"));
    AliWarning(Form("  freeze-out criteria Tf=145 MeV"));
    AliWarning(Form("  for details see $(ALICE_ROOT)/TTherminator/data/LHC500C0005/FO.txt"));
    fFOHSlocation.Append(aroot);
    fFOHSlocation.Append("/TTherminator/data/LHC500C0005");
  }
  if (strstr(set, "LHC500C0510")) {
    AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
    AliWarning(Form("  Pb-Pb collisions, centrality 5-10 percent"));
    AliWarning(Form("  initial temperature at tau=1 fm in the center Ti=500 MeV"));
    AliWarning(Form("  freeze-out criteria Tf=145 MeV"));
    AliWarning(Form("  for details see $(ALICE_ROOT)/TTherminator/data/LHC500C0510/FO.txt"));
    fFOHSlocation.Append(aroot);
    fFOHSlocation.Append("/TTherminator/data/LHC500C0510");
  }
  if (strstr(set, "LHC500C1020")) {
    AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
    AliWarning(Form("  Pb-Pb collisions, centrality 10-20 percent"));
    AliWarning(Form("  initial temperature at tau=1 fm in the center Ti=500 MeV"));
    AliWarning(Form("  freeze-out criteria Tf=145 MeV"));
    AliWarning(Form("  for details see $(ALICE_ROOT)/TTherminator/data/LHC500C1020/FO.txt"));
    fFOHSlocation.Append(aroot);
    fFOHSlocation.Append("/TTherminator/data/LHC500C1020");
  }
  else if (strstr(set, "LHC500C2030")) {
    AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
    AliWarning(Form("  Pb-Pb collisions, centrality 20-30 percent"));
    AliWarning(Form("  initial temperature at tau=1 fm in the center Ti=500 MeV"));
    AliWarning(Form("  freeze-out criteria Tf=145 MeV"));
    AliWarning(Form("  for details see $(ALICE_ROOT)/TTherminator/data/LHC500C2030/FO.txt"));
    fFOHSlocation.Append(aroot);
    fFOHSlocation.Append("/TTherminator/data/LHC500C2030");
  }
  else if (strstr(set, "LHC500C3040")) {
    AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
    AliWarning(Form("  Pb-Pb collisions, centrality 30-40 percent"));
    AliWarning(Form("  initial temperature at tau=1 fm in the center Ti=500 MeV"));
    AliWarning(Form("  freeze-out criteria Tf=145 MeV"));
    AliWarning(Form("  for details see $(ALICE_ROOT)/TTherminator/data/LHC500C3040/FO.txt"));
    fFOHSlocation.Append(aroot);
    fFOHSlocation.Append("/TTherminator/data/LHC500C3040");
  }
  else if (strstr(set, "LHC500C4050")) {
    AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
    AliWarning(Form("  Pb-Pb collisions, centrality 40-50 percent"));
    AliWarning(Form("  initial temperature at tau=1 fm in the center Ti=500 MeV"));
    AliWarning(Form("  freeze-out criteria Tf=145 MeV"));
    AliWarning(Form("  for details see $(ALICE_ROOT)/TTherminator/data/LHC500C4050/FO.txt"));
    fFOHSlocation.Append(aroot);
    fFOHSlocation.Append("/TTherminator/data/LHC500C4050");
  }
  else if (strstr(set, "LHC276TC0005")) {
    AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
    AliWarning(Form("  Pb-Pb collisions, 2.76TeV, centrality 00-05 percent"));
    AliWarning(Form("  freeze-out criteria Tf=145 MeV"));
    AliWarning(Form("  for details see $(ALICE_ROOT)/TTherminator/data/LHC276TC0005/FO.txt"));
    fFOHSlocation.Append(aroot);
    fFOHSlocation.Append("/TTherminator/data/LHC276TC0005");
  }
  else if (strstr(set, "LHC276TC0510")) {
    AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
    AliWarning(Form("  Pb-Pb collisions, 2.76TeV, centrality 05-10 percent"));
    AliWarning(Form("  freeze-out criteria Tf=145 MeV"));
    AliWarning(Form("  for details see $(ALICE_ROOT)/TTherminator/data/LHC276TC0510/FO.txt"));
    fFOHSlocation.Append(aroot);
    fFOHSlocation.Append("/TTherminator/data/LHC276TC0510");
  }
  else if (strstr(set, "LHC276TC1020")) {
    AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
    AliWarning(Form("  Pb-Pb collisions, 2.76TeV, centrality 10-20 percent"));
    AliWarning(Form("  freeze-out criteria Tf=145 MeV"));
    AliWarning(Form("  for details see $(ALICE_ROOT)/TTherminator/data/LHC276TC1020/FO.txt"));
    fFOHSlocation.Append(aroot);
    fFOHSlocation.Append("/TTherminator/data/LHC276TC1020");
  }
  else if (strstr(set, "LHC276TC2030")) {
    AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
    AliWarning(Form("  Pb-Pb collisions, 2.76TeV, centrality 20-30 percent"));
    AliWarning(Form("  freeze-out criteria Tf=145 MeV"));
    AliWarning(Form("  for details see $(ALICE_ROOT)/TTherminator/data/LHC276TC2030/FO.txt"));
    fFOHSlocation.Append(aroot);
    fFOHSlocation.Append("/TTherminator/data/LHC276TC2030");
  }
  else if (strstr(set, "LHC276TC3040")) {
    AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
    AliWarning(Form("  Pb-Pb collisions, 2.76TeV, centrality 30-40 percent"));
    AliWarning(Form("  freeze-out criteria Tf=145 MeV"));
    AliWarning(Form("  for details see $(ALICE_ROOT)/TTherminator/data/LHC276TC3040/FO.txt"));
    fFOHSlocation.Append(aroot);
    fFOHSlocation.Append("/TTherminator/data/LHC276TC3040");
  }
  else if (strstr(set, "LHC276TC4050")) {
    AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
    AliWarning(Form("  Pb-Pb collisions, 2.76TeV, centrality 40-50 percent"));
    AliWarning(Form("  freeze-out criteria Tf=145 MeV"));
    AliWarning(Form("  for details see $(ALICE_ROOT)/TTherminator/data/LHC276TC4050/FO.txt"));
    fFOHSlocation.Append(aroot);
    fFOHSlocation.Append("/TTherminator/data/LHC276TC4050");
  }
  else {
    AliWarning(Form("Did not find Lhyquid set %s", set));
    AliWarning(Form("Reverting to default: current directory"));
    fFOHSlocation += "";
  }

}

void AliGenTherminator::SetLhyquidInputDir(const char *inputdir)
{
  // Select Your own Lhyquid hypersurface
  fFOHSlocation = inputdir;
}
Commit	Line	Data
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
ae89e57a	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
	24	ClassImp(AliGenTherminator)
	25
	26	using namespace std;
	27
	28	AliGenTherminator::AliGenTherminator():
	29	AliGenMC(),
	30	fNt(0),
	31	fEventNumber(0),
	32	fFileName(""),
	33	fFreezeOutModel(""),
9368eb75	34	fFOHSlocation(""),
2e967919	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
9368eb75	48	fFOHSlocation = "";
	49
	50	fEnergyCMS = 5500.;
	51	fAProjectile = 208;
	52	fZProjectile = 82;
	53	fProjectile = "A";
	54	fATarget = 208;
	55	fZTarget = 82;
	56	fTarget = "A";
2e967919	57	}
	58	AliGenTherminator::AliGenTherminator(Int_t npart):
	59	AliGenMC(npart),
	60	fNt(0),
	61	fEventNumber(0),
	62	fFileName(""),
	63	fFreezeOutModel(""),
9368eb75	64	fFOHSlocation(""),
2e967919	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
2e967919	78	fNprimaries = 0;
9368eb75	79	fEnergyCMS = 5500.;
	80	fAProjectile = 208;
	81	fZProjectile = 82;
	82	fProjectile = "A";
	83	fATarget = 208;
	84	fZTarget = 82;
	85	fTarget = "A";
2e967919	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};
21391258	101	Float_t time = 0.;
2e967919	102	Float_t origin0[3] = {0,0,0};
21391258	103	Float_t time0 = 0.;
2e967919	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];
21391258	113	time0 = fTime;
2e967919	114
	115	// Generate one event
	116
	117	((TTherminator *) fMCEvGen)->GenerateEvent();
	118	AliWarning("Generated");
8507138f	119	((TTherminator *) fMCEvGen)->ImportParticles(&fParticles);
2e967919	120
8507138f	121	Int_t np = fParticles.GetEntriesFast();
2e967919	122	AliWarning(Form("Imported %d particles", np));
2e967919	123
deb008f3	124	Int_t *idsOnStack;
	125	idsOnStack = new Int_t[np];
	126
2e967919	127	TParticle *iparticle;
9abb118d	128	Double_t evrot = gRandom->Rndm()*TMath::Pi();
2e967919	129
2e967919	130	for (int i = 0; i < np; i++) {
8507138f	131	iparticle = (TParticle *) fParticles.At(i);
2e967919	132	Bool_t hasMother = (iparticle->GetFirstMother() >=0);
	133	Bool_t hasDaughter = (iparticle->GetFirstDaughter() >=0);
	134
c9a371fc	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(massmass + 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());
c9c8970c	155	origin[0] = vrho*TMath::Cos(vphi + evrot);
	156	origin[1] = vrho*TMath::Sin(vphi + evrot);
	157	origin[2] = iparticle->Vz();
21391258	158	time = iparticle->T();
21391258	159
c9a371fc	160	imo = -1;
4e626b54	161	// TParticle* mother = 0;
c9a371fc	162	if (hasMother) {
c9a371fc	163	imo = iparticle->GetFirstMother();
20c79170	164	// mother = (TParticle *) fParticles.At(imo);
c9a371fc	165	} // if has mother
	166	Bool_t tFlag = (!hasDaughter);
	167
9368eb75	168	// printf("Pushing Track %d with status %d mother %d\n", kf, tFlag, imo>=0?idsOnStack[imo]:imo);
d75811a6	169	PushTrack(tFlag,imo>=0?idsOnStack[imo]:imo,kf,
c9a371fc	170	p[0],p[1],p[2],energy,
21391258	171	origin[0],origin[1],origin[2],time,
c9a371fc	172	polar[0],polar[1],polar[2],
c9a371fc	173	hasMother ? kPDecay:kPNoProcess,nt);
deb008f3	174	idsOnStack[i] = nt;
c9a371fc	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(massmass + 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();
21391258	204	time = iparticle->T();
21391258	205
c9a371fc	206	imo = -1;
6750cf94	207	// TParticle* mother = 0;
c9a371fc	208	if (hasMother) {
c9a371fc	209	imo = iparticle->GetFirstMother();
6750cf94	210	// mother = (TParticle *) fParticles.At(imo);
c9a371fc	211	} // if has mother
	212	Bool_t tFlag = (hasDaughter);
	213
9368eb75	214	// printf("Found mother %i with true id %i\n", imo, imo>=0?idsOnStack[imo]:imo);
9368eb75	215	// printf("Pushing Track %d with status %d mother %d\n", kf, tFlag, imo>=0?idsOnStack[imo]:imo);
d75811a6	216	PushTrack(tFlag,imo>=0?idsOnStack[imo]:imo,kf,
c9a371fc	217	p[0],p[1],p[2],energy,
21391258	218	origin[0],origin[1],origin[2],time,
c9a371fc	219	polar[0],polar[1],polar[2],
c9a371fc	220	hasMother ? kPDecay:kPNoProcess,nt);
deb008f3	221	idsOnStack[i] = nt;
c9a371fc	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();
21391258	228	time = time0+iparticle->T();
c9a371fc	229
c9a371fc	230	imo = nt;
deb008f3	231	// mother = (TParticle *) fParticles.At(nt);
c9a371fc	232	tFlag = (!hasDaughter);
c9a371fc	233
9368eb75	234	// printf("Pushing Track %d with status %d mother %d\n", kf, tFlag, imo);
c9a371fc	235	PushTrack(tFlag,imo,kf,
c9a371fc	236	p[0],p[1],p[2],energy,
21391258	237	origin[0],origin[1],origin[2],time,
c9a371fc	238	polar[0],polar[1],polar[2],
	239	hasMother ? kPDecay:kPNoProcess,nt);
	240	fNprimaries++;
	241	KeepTrack(nt);
	242	}
2e967919	243	}
2e967919	244
deb008f3	245
deb008f3	246
2e967919	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];
21391258	254	Float_t eventTime = time0;
2e967919	255
9abb118d	256	// Builds the event header, to be called after each event
	257	AliGenEventHeader* header = new AliGenHijingEventHeader("Therminator");
	258
2e967919	259	// Header
9abb118d	260	// AliGenEventHeader* header = new AliGenEventHeader("Therminator");
2e967919	261	// Event Vertex
9abb118d	262	// header->SetPrimaryVertex(eventVertex);
	263	// header->SetNProduced(fNprimaries);
	264
	265	((AliGenHijingEventHeader*) header)->SetNProduced(fNprimaries);
	266	((AliGenHijingEventHeader*) header)->SetPrimaryVertex(eventVertex);
21391258	267	((AliGenHijingEventHeader*) header)->SetInteractionTime(eventTime);
9abb118d	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);
21391258	304	header->SetInteractionTime(fTime);
9abb118d	305	AddHeader(header);
	306	fCollisionGeometry = (AliGenHijingEventHeader*) header;
	307
4ce766eb	308	delete [] idsOnStack;
deb008f3	309
c9a371fc	310	// gAlice->SetGenEventHeader(header);
2e967919	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"));
ae89e57a	349
	350	TString aroot = gSystem->Getenv("ALICE_ROOT");
	351	ifstream in((aroot+"/TTherminator/data/SHARE/particles.data").Data());
	352	// ifstream in("particles.data");
2e967919	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
ae89e57a	399	const char *aroot = gSystem->Getenv("ALICE_ROOT");
2e967919	400	ofstream *ostr = new ofstream("therminator.in");
	401	(*ostr) << "NumberOfEvents = 1" << endl;
	402	(*ostr) << "Randomize = 1" << endl;
	403	(*ostr) << "TableType = SHARE" << endl;
ae89e57a	404	(*ostr) << "InputDirSHARE = "<< aroot << "/TTherminator/data/SHARE" << endl;
2e967919	405	(*ostr) << "EventOutputFile = " << fFileName.Data() << endl;
ae89e57a	406	(*ostr) << "FOHSLocation = " << fFOHSlocation.Data() << endl;
2e967919	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;
9368eb75	424	AliWarning(Form("Selected model %s", fFreezeOutModel.Data()));
9368eb75	425	AliWarning(Form("FOHSLocation is %s", fFOHSlocation.Data()));
2e967919	426	}
ae89e57a	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"));
7d4c52b2	434	AliWarning(Form(" Pb-Pb collisions, centrality 0-5 percent"));
ae89e57a	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"));
9368eb75	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");
ae89e57a	458	}
	459	else if (strstr(set, "LHC500C2030")) {
	460	AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
7d4c52b2	461	AliWarning(Form(" Pb-Pb collisions, centrality 20-30 percent"));
ae89e57a	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"));
9368eb75	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");
ae89e57a	485	}
269312f5	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	}
ae89e57a	534	else {
	535	AliWarning(Form("Did not find Lhyquid set %s", set));
	536	AliWarning(Form("Reverting to default: current directory"));
9368eb75	537	fFOHSlocation += "";
ae89e57a	538	}
9368eb75	539
ae89e57a	540	}
	541
	542	void AliGenTherminator::SetLhyquidInputDir(const char *inputdir)
	543	{
	544	// Select Your own Lhyquid hypersurface
	545	fFOHSlocation = inputdir;
	546	}