[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(""),
  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
  fParticles = new TClonesArray("TParticle",20000);    
}
AliGenTherminator::AliGenTherminator(Int_t npart):
  AliGenMC(npart),
  fNt(0),
  fEventNumber(0),
  fFileName(""),
  fFreezeOutModel(""),
  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
  fParticles = new TClonesArray("TParticle",20000);    
  fNprimaries = 0;
}

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 origin0[3]  =   {0,0,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];

  // Generate one event

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

  Int_t np = fParticles->GetEntriesFast();
  AliWarning(Form("Imported %d particles", 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);
    
    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] = origin0[0]+vrho*TMath::Cos(vphi + evrot);
    origin[1] = origin0[1]+vrho*TMath::Sin(vphi + evrot);
    origin[2] = origin0[2]+iparticle->Vz();
	      
    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);
    PushTrack(tFlag,imo,kf,
	      p[0],p[1],p[2],energy,
	      origin[0],origin[1],origin[2],iparticle->T()*0.197327*1e-13/300000000,
	      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];

// 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)->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);
  AddHeader(header);
  fCollisionGeometry = (AliGenHijingEventHeader*)  header;

  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;
}

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%"));
    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 = aroot;
    fFOHSlocation += "/TTherminator/data/LHC500C0005";
  }
  else if (strstr(set, "LHC500C2030")) {
    AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
    AliWarning(Form("  Pb-Pb collisions, centrality 20-30%"));
    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 = aroot;
    fFOHSlocation += "/TTherminator/data/LHC500C2030";
  }
  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(""),
	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;
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(massmass + p[0]p[0] + p[1]p[1] + p[2]p[2]);
e2867a14	124
	125	Double_t vphi = TMath::ATan2(iparticle->Vy(), iparticle->Vx());
	126	Double_t vrho = TMath::Hypot(iparticle->Vx(), iparticle->Vy());
	127	origin[0] = origin0[0]+vrho*TMath::Cos(vphi + evrot);
	128	origin[1] = origin0[1]+vrho*TMath::Sin(vphi + evrot);
2e967919	129	origin[2] = origin0[2]+iparticle->Vz();
	130
	131	imo = -1;
	132	TParticle* mother = 0;
	133	if (hasMother) {
	134	imo = iparticle->GetFirstMother();
	135	mother = (TParticle *) fParticles->At(imo);
	136	} // if has mother
	137	Bool_t tFlag = (!hasDaughter);
	138
	139	printf("Pushing Track %d with status %d mother %d\n", kf, tFlag, imo);
	140	PushTrack(tFlag,imo,kf,
	141	p[0],p[1],p[2],energy,
	142	origin[0],origin[1],origin[2],iparticle->T()0.1973271e-13/300000000,
	143	polar[0],polar[1],polar[2],
	144	hasMother ? kPDecay:kPNoProcess,nt);
	145
	146	fNprimaries++;
	147	KeepTrack(nt);
	148	}
	149
	150	SetHighWaterMark(fNprimaries);
	151
	152	TArrayF eventVertex;
	153	eventVertex.Set(3);
	154	eventVertex[0] = origin0[0];
	155	eventVertex[1] = origin0[1];
	156	eventVertex[2] = origin0[2];
	157
9abb118d	158	// Builds the event header, to be called after each event
	159	AliGenEventHeader* header = new AliGenHijingEventHeader("Therminator");
	160
2e967919	161	// Header
9abb118d	162	// AliGenEventHeader* header = new AliGenEventHeader("Therminator");
2e967919	163	// Event Vertex
9abb118d	164	// header->SetPrimaryVertex(eventVertex);
	165	// header->SetNProduced(fNprimaries);
	166
	167	((AliGenHijingEventHeader*) header)->SetNProduced(fNprimaries);
	168	((AliGenHijingEventHeader*) header)->SetPrimaryVertex(eventVertex);
	169	((AliGenHijingEventHeader*) header)->SetImpactParameter(0.0);
	170	((AliGenHijingEventHeader*) header)->SetTotalEnergy(0.0);
	171	((AliGenHijingEventHeader*) header)->SetHardScatters(0);
	172	((AliGenHijingEventHeader*) header)->SetParticipants(0, 0);
	173	((AliGenHijingEventHeader*) header)->SetCollisions(0, 0, 0, 0);
	174	((AliGenHijingEventHeader*) header)->SetSpectators(0, 0, 0, 0);
	175	((AliGenHijingEventHeader*) header)->SetReactionPlaneAngle(evrot);
	176
	177
	178	// 4-momentum vectors of the triggered jets.
	179	//
	180	// Before final state gluon radiation.
	181	// TLorentzVector* jet1 = new TLorentzVector(fHijing->GetHINT1(21),
	182	// fHijing->GetHINT1(22),
	183	// fHijing->GetHINT1(23),
	184	// fHijing->GetHINT1(24));
	185
	186	// TLorentzVector* jet2 = new TLorentzVector(fHijing->GetHINT1(31),
	187	// fHijing->GetHINT1(32),
	188	// fHijing->GetHINT1(33),
	189	// fHijing->GetHINT1(34));
	190	// // After final state gluon radiation.
	191	// TLorentzVector* jet3 = new TLorentzVector(fHijing->GetHINT1(26),
	192	// fHijing->GetHINT1(27),
	193	// fHijing->GetHINT1(28),
	194	// fHijing->GetHINT1(29));
	195
	196	// TLorentzVector* jet4 = new TLorentzVector(fHijing->GetHINT1(36),
	197	// fHijing->GetHINT1(37),
	198	// fHijing->GetHINT1(38),
	199	// fHijing->GetHINT1(39));
	200	// ((AliGenHijingEventHeader*) header)->SetJets(jet1, jet2, jet3, jet4);
	201	// Bookkeeping for kinematic bias
	202	// ((AliGenHijingEventHeader*) header)->SetTrials(fTrials);
	203	// Event Vertex
	204	header->SetPrimaryVertex(fVertex);
	205	AddHeader(header);
	206	fCollisionGeometry = (AliGenHijingEventHeader*) header;
	207
2e967919	208	gAlice->SetGenEventHeader(header);
	209	}
	210
	211	void AliGenTherminator::Init()
	212	{
	213	// Initialize global variables and
	214	// particle and decay tables
	215	if (fFileName.Length() == 0)
	216	fFileName = "event.out";
	217	ReadShareParticleTable();
	218
	219	SetMC(new TTherminator());
	220
	221	AliGenMC::Init();
	222	((TTherminator *) fMCEvGen)->Initialize();
	223	}
	224	void AliGenTherminator::SetFileName(const char *infilename)
	225	{
	226	// Set parameter filename
	227	fFileName = infilename;
	228	}
	229	void AliGenTherminator::SetEventNumberInFile(int evnum)
	230	{
	231	// Set number of events to generate - default: 1
	232	fEventNumber = evnum;
	233	}
	234
	235	void AliGenTherminator::ReadShareParticleTable()
	236	{
	237	// Read in particle table from share
	238	// and add missing particle type to TDatabasePDG
	239
	240	char str[50];
	241	char str1[200];
	242
	243	TDatabasePDG *tInstance = TDatabasePDG::Instance();
	244	TParticlePDG *tParticleType;
	245
	246	AliWarning(Form("Reading particle types from particles.data"));
ae89e57a	247
	248	TString aroot = gSystem->Getenv("ALICE_ROOT");
	249	ifstream in((aroot+"/TTherminator/data/SHARE/particles.data").Data());
	250	// ifstream in("particles.data");
2e967919	251
	252	int charge;
	253
	254	int number=0;
	255	if ((in) && (in.is_open()))
	256	{
	257	//START OF HEAD-LINE
	258	in.ignore(200,'\n');
	259	in.ignore(200,'\n');
	260	in.ignore(200,'\n');
	261	//END OF HEAD-LINE
	262
	263	while (in>>str)
	264	{
	265	if (/(str == '#')\|\|/(str<65)\|\|(*str>122))
	266	{
	267	in.getline(str1,200);
	268	continue;
	269	}
	270	double mass, gamma, spin, tI3, tI, q, s, aq, as, c, ac, mc;
	271
	272	in>>mass>>gamma>>spin>>tI>>tI3>>q>>s>>aq>>as>>c>>ac>>mc;
	273	number++;
	274	tParticleType = tInstance->GetParticle((int) mc);
	275	if (!tParticleType) {
	276	charge = 0;
	277	if (strstr(str, "plu")) charge = 1;
	278	if (strstr(str, "min")) charge = -1;
	279	if (strstr(str, "plb")) charge = -1;
	280	if (strstr(str, "mnb")) charge = 1;
	281	if (strstr(str, "plp")) charge = 2;
	282	if (strstr(str, "ppb")) charge = -2;
	283	tInstance->AddParticle(str, str, mass, gamma == 0.0 ? 1:0, gamma, charge , "meson", (int) mc);
	284	AliWarning(Form("Added particle %s with PDG PID %d charge %d", str, (int) mc, charge));
	285	// 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));
	286
	287	}
	288	}
	289	in.close();
	290	}
	291	CreateTherminatorInputFile();
	292	}
	293
	294	void AliGenTherminator::CreateTherminatorInputFile()
	295	{
	296	// Create Therminator input file
ae89e57a	297	const char *aroot = gSystem->Getenv("ALICE_ROOT");
2e967919	298	ofstream *ostr = new ofstream("therminator.in");
	299	(*ostr) << "NumberOfEvents = 1" << endl;
	300	(*ostr) << "Randomize = 1" << endl;
	301	(*ostr) << "TableType = SHARE" << endl;
ae89e57a	302	(*ostr) << "InputDirSHARE = "<< aroot << "/TTherminator/data/SHARE" << endl;
2e967919	303	(*ostr) << "EventOutputFile = " << fFileName.Data() << endl;
ae89e57a	304	(*ostr) << "FOHSLocation = " << fFOHSlocation.Data() << endl;
2e967919	305	(*ostr) << "FreezeOutModel = " << fFreezeOutModel.Data() << endl;
	306	(*ostr) << "BWVt = " << fBWVt << endl;
	307	(*ostr) << "Tau = " << fTau << endl;
	308	(*ostr) << "RhoMax = " << fRhoMax << endl;
	309	(*ostr) << "Temperature = " << fTemperature << endl;
	310	(*ostr) << "MiuI = " << fMiuI << endl;
	311	(*ostr) << "MiuS = " << fMiuS << endl;
	312	(*ostr) << "MiuB = " << fMiuB << endl;
	313	(*ostr) << "AlphaRange = " << fAlfaRange << endl;
	314	(*ostr) << "RapidityRange = " << fRapRange << endl;
	315	(*ostr) << "NumberOfIntegrateSamples = 1000000" << endl;
	316	}
	317
	318	void AliGenTherminator::SetModel(const char *model)
	319	{
	320	// Set the freeze-out model to use
	321	fFreezeOutModel = model;
	322	}
ae89e57a	323
	324	void AliGenTherminator::SetLhyquidSet(const char *set)
	325	{
	326	// Select one of pregenerated Lhyquid hypersurfaces
	327	const char *aroot = gSystem->Getenv("ALICE_ROOT");
	328	if (strstr(set, "LHC500C0005")) {
	329	AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
	330	AliWarning(Form(" Pb-Pb collisions, centrality 0-5%"));
	331	AliWarning(Form(" initial temperature at tau=1 fm in the center Ti=500 MeV"));
	332	AliWarning(Form(" freeze-out criteria Tf=145 MeV"));
	333	AliWarning(Form(" for details see $(ALICE_ROOT)/TTherminator/data/LHC500C0005/FO.txt"));
	334	fFOHSlocation = aroot;
	335	fFOHSlocation += "/TTherminator/data/LHC500C0005";
	336	}
	337	else if (strstr(set, "LHC500C2030")) {
	338	AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
	339	AliWarning(Form(" Pb-Pb collisions, centrality 20-30%"));
	340	AliWarning(Form(" initial temperature at tau=1 fm in the center Ti=500 MeV"));
	341	AliWarning(Form(" freeze-out criteria Tf=145 MeV"));
	342	AliWarning(Form(" for details see $(ALICE_ROOT)/TTherminator/data/LHC500C2030/FO.txt"));
	343	fFOHSlocation = aroot;
	344	fFOHSlocation += "/TTherminator/data/LHC500C2030";
	345	}
	346	else {
	347	AliWarning(Form("Did not find Lhyquid set %s", set));
	348	AliWarning(Form("Reverting to default: current directory"));
	349	fFOHSlocation = "";
	350	}
	351	}
	352
	353	void AliGenTherminator::SetLhyquidInputDir(const char *inputdir)
	354	{
	355	// Select Your own Lhyquid hypersurface
	356	fFOHSlocation = inputdir;
	357	}