[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
// 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 "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;
  
  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();
    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]);
    origin[0] = origin0[0]+iparticle->Vx();
    origin[1] = origin0[1]+iparticle->Vy();
    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];

  // Header
  AliGenEventHeader* header = new AliGenEventHeader("Therminator");
  // Event Vertex
  header->SetPrimaryVertex(eventVertex);
  header->SetNProduced(fNprimaries);
  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"));
  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
  ofstream *ostr = new ofstream("therminator.in");
  (*ostr) << "NumberOfEvents = 1" << endl;
  (*ostr) << "Randomize = 1" << endl;
  (*ostr) << "TableType = SHARE" << endl;
  (*ostr) << "InputDirSHARE = ." << endl;
  (*ostr) << "EventOutputFile = " << fFileName.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;
}
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
	4	// Author: Adam.Kisiel@cern.ch
	5
	6	#include <iostream>
	7	#include <fstream>
	8	#include <sstream>
	9	#include <TClonesArray.h>
	10	#include <TMCProcess.h>
	11	#include <TDatabasePDG.h>
	12	#include <TParticle.h>
	13
	14	#include "AliConst.h"
	15	#include "AliDecayer.h"
	16	#include "AliGenEventHeader.h"
	17	#include "AliGenTherminator.h"
	18	#include "AliLog.h"
	19	#include "AliRun.h"
	20
	21	ClassImp(AliGenTherminator)
	22
	23	using namespace std;
	24
	25	AliGenTherminator::AliGenTherminator():
	26	AliGenMC(),
	27	fNt(0),
	28	fEventNumber(0),
	29	fFileName(""),
	30	fFreezeOutModel(""),
	31	fTemperature(0.1656),
	32	fMiuI(-0.0009),
	33	fMiuS(0.0),
	34	fMiuB(0.0008),
	35	fAlfaRange(6.0),
	36	fRapRange(4.0),
	37	fRhoMax(8.0),
	38	fTau(8.0),
	39	fBWA(0.0),
	40	fBWVt(1.41),
	41	fBWDelay(0.0)
	42	{
	43	// Default constructor
	44	fParticles = new TClonesArray("TParticle",20000);
	45	}
	46	AliGenTherminator::AliGenTherminator(Int_t npart):
	47	AliGenMC(npart),
	48	fNt(0),
	49	fEventNumber(0),
	50	fFileName(""),
	51	fFreezeOutModel(""),
	52	fTemperature(0.1656),
	53	fMiuI(-0.0009),
	54	fMiuS(0.0),
	55	fMiuB(0.0008),
	56	fAlfaRange(6.0),
	57	fRapRange(4.0),
	58	fRhoMax(8.0),
	59	fTau(8.0),
	60	fBWA(0.0),
	61	fBWVt(1.41),
	62	fBWDelay(0.0)
	63	{
	64	// Constructor specifying the size of the particle table
65	fParticles = new TClonesArray("TParticle",20000);
66	fNprimaries = 0;
67	}
68
69	AliGenTherminator::~AliGenTherminator()
70	{
71	// AliGenMC::~AliGenMC();
72	// if (fTherminator) delete fTherminator;
73	}
74
75	void AliGenTherminator::Generate()
76	{
77	// Run single event generation with the Therminator model
78	AliWarning("Generating event from AliGenTherminator");
79
80	Float_t polar[3] = {0,0,0};
81	Float_t origin[3] = {0,0,0};
82	Float_t origin0[3] = {0,0,0};
83	Float_t p[3];
84	Float_t mass, energy;
85
86	Int_t nt = 0;
87	Int_t j, kf, ks, imo;
88	kf = 0;
89
90	Vertex();
91	for (j=0; j < 3; j++) origin0[j] = fVertex[j];
92
93	// Generate one event
94
95	((TTherminator *) fMCEvGen)->GenerateEvent();
96	AliWarning("Generated");
97	((TTherminator *) fMCEvGen)->ImportParticles(fParticles);
98
99	Int_t np = fParticles->GetEntriesFast();
100	AliWarning(Form("Imported %d particles", np));
101
102	TParticle *iparticle;
103
104	for (int i = 0; i < np; i++) {
105	iparticle = (TParticle *) fParticles->At(i);
106	Bool_t hasMother = (iparticle->GetFirstMother() >=0);
107	Bool_t hasDaughter = (iparticle->GetFirstDaughter() >=0);
108
109	kf = iparticle->GetPdgCode();
110	ks = iparticle->GetStatusCode();
111	p[0] = iparticle->Px();
112	p[1] = iparticle->Py();
113	p[2] = iparticle->Pz();
114	mass = TDatabasePDG::Instance()->GetParticle(kf)->Mass();
115	energy = sqrt(massmass + p[0]p[0] + p[1]p[1] + p[2]p[2]);
116	origin[0] = origin0[0]+iparticle->Vx();
117	origin[1] = origin0[1]+iparticle->Vy();
118	origin[2] = origin0[2]+iparticle->Vz();
119
120	imo = -1;
121	TParticle* mother = 0;
122	if (hasMother) {
123	imo = iparticle->GetFirstMother();
124	mother = (TParticle *) fParticles->At(imo);
125	} // if has mother
126	Bool_t tFlag = (!hasDaughter);
127
128	printf("Pushing Track %d with status %d mother %d\n", kf, tFlag, imo);
129	PushTrack(tFlag,imo,kf,
130	p[0],p[1],p[2],energy,
131	origin[0],origin[1],origin[2],iparticle->T()0.1973271e-13/300000000,
132	polar[0],polar[1],polar[2],
133	hasMother ? kPDecay:kPNoProcess,nt);
134
135	fNprimaries++;
136	KeepTrack(nt);
137	}
138
139	SetHighWaterMark(fNprimaries);
140
141	TArrayF eventVertex;
142	eventVertex.Set(3);
143	eventVertex[0] = origin0[0];
144	eventVertex[1] = origin0[1];
145	eventVertex[2] = origin0[2];
146
147	// Header
148	AliGenEventHeader* header = new AliGenEventHeader("Therminator");
149	// Event Vertex
150	header->SetPrimaryVertex(eventVertex);
151	header->SetNProduced(fNprimaries);
152	gAlice->SetGenEventHeader(header);
153	}
154
155	void AliGenTherminator::Init()
156	{
157	// Initialize global variables and
158	// particle and decay tables
159	if (fFileName.Length() == 0)
160	fFileName = "event.out";
161	ReadShareParticleTable();
162
163	SetMC(new TTherminator());
164
165	AliGenMC::Init();
166	((TTherminator *) fMCEvGen)->Initialize();
167	}
168	void AliGenTherminator::SetFileName(const char *infilename)
169	{
170	// Set parameter filename
171	fFileName = infilename;
172	}
173	void AliGenTherminator::SetEventNumberInFile(int evnum)
174	{
175	// Set number of events to generate - default: 1
176	fEventNumber = evnum;
177	}
178
179	void AliGenTherminator::ReadShareParticleTable()
180	{
181	// Read in particle table from share
182	// and add missing particle type to TDatabasePDG
183
184	char str[50];
185	char str1[200];
186
187	TDatabasePDG *tInstance = TDatabasePDG::Instance();
188	TParticlePDG *tParticleType;
189
190	AliWarning(Form("Reading particle types from particles.data"));
191	ifstream in("particles.data");
192
193	int charge;
194
195	int number=0;
196	if ((in) && (in.is_open()))
197	{
198	//START OF HEAD-LINE
199	in.ignore(200,'\n');
200	in.ignore(200,'\n');
201	in.ignore(200,'\n');
202	//END OF HEAD-LINE
203
204	while (in>>str)
205	{
206	if (/(str == '#')\|\|/(str<65)\|\|(*str>122))
207	{
208	in.getline(str1,200);
209	continue;
210	}
211	double mass, gamma, spin, tI3, tI, q, s, aq, as, c, ac, mc;
212
213	in>>mass>>gamma>>spin>>tI>>tI3>>q>>s>>aq>>as>>c>>ac>>mc;
214	number++;
215	tParticleType = tInstance->GetParticle((int) mc);
216	if (!tParticleType) {
217	charge = 0;
218	if (strstr(str, "plu")) charge = 1;
219	if (strstr(str, "min")) charge = -1;
220	if (strstr(str, "plb")) charge = -1;
221	if (strstr(str, "mnb")) charge = 1;
222	if (strstr(str, "plp")) charge = 2;
223	if (strstr(str, "ppb")) charge = -2;
224	tInstance->AddParticle(str, str, mass, gamma == 0.0 ? 1:0, gamma, charge , "meson", (int) mc);
225	AliWarning(Form("Added particle %s with PDG PID %d charge %d", str, (int) mc, charge));
226	// 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));
227
228	}
229	}
230	in.close();
231	}
232	CreateTherminatorInputFile();
233	}
234
235	void AliGenTherminator::CreateTherminatorInputFile()
236	{
237	// Create Therminator input file
238	ofstream *ostr = new ofstream("therminator.in");
239	(*ostr) << "NumberOfEvents = 1" << endl;
240	(*ostr) << "Randomize = 1" << endl;
241	(*ostr) << "TableType = SHARE" << endl;
242	(*ostr) << "InputDirSHARE = ." << endl;
243	(*ostr) << "EventOutputFile = " << fFileName.Data() << endl;
244	(*ostr) << "FreezeOutModel = " << fFreezeOutModel.Data() << endl;
245	(*ostr) << "BWVt = " << fBWVt << endl;
246	(*ostr) << "Tau = " << fTau << endl;
247	(*ostr) << "RhoMax = " << fRhoMax << endl;
248	(*ostr) << "Temperature = " << fTemperature << endl;
249	(*ostr) << "MiuI = " << fMiuI << endl;
250	(*ostr) << "MiuS = " << fMiuS << endl;
251	(*ostr) << "MiuB = " << fMiuB << endl;
252	(*ostr) << "AlphaRange = " << fAlfaRange << endl;
253	(*ostr) << "RapidityRange = " << fRapRange << endl;
254	(*ostr) << "NumberOfIntegrateSamples = 1000000" << endl;
255	}
256
257	void AliGenTherminator::SetModel(const char *model)
258	{
259	// Set the freeze-out model to use
260	fFreezeOutModel = model;
261	}