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

///////////////////////////////////////////////////////////////////////////////
//                                                                           //
// TTherminator: global interface class to the Therminator model.            //
// Initialized global variables and runs the event generation                //
//                                                                           //
///////////////////////////////////////////////////////////////////////////////
/******************************************************************************
 *                      T H E R M I N A T O R                                 *
 *                   THERMal heavy-IoN generATOR                              *
 *                           version 1.0                                      *
 *                                                                            *
 * Authors of the model: Wojciech Broniowski, Wojciech.Broniowski@ifj.edu.pl, *
 *                       Wojciech Florkowski, Wojciech.Florkowski@ifj.edu.pl  *
 * Authors of the code:  Adam Kisiel, kisiel@if.pw.edu.pl                     *
 *                       Tomasz Taluc, ttaluc@if.pw.edu.pl                    *
 * Code designers: Adam Kisiel, Tomasz Taluc, Wojciech Broniowski,            *
 *                 Wojciech Florkowski                                        *
 *                                                                            *
 * For the detailed description of the program and furhter references         * 
 * to the description of the model plesase refer to: nucl-th/0504047,         *
 * accessibile at: http://www.arxiv.org/nucl-th/0504047                       *
 *                                                                            *
 * Homepage: http://hirg.if.pw.edu.pl/en/therminator/                         *
 *                                                                            *
 * This code can be freely used and redistributed. However if you decide to   *
 * make modifications to the code, please contact the authors, especially     *
 * if you plan to publish the results obtained with such modified code.       *
 * Any publication of results obtained using this code must include the       *
 * reference to nucl-th/0504047 and the published version of it, when         *
 * available.                                                                 *
 *                                                                            *
 *****************************************************************************/
#include <TTherminator.h>
#include <fstream>
#include <TDatabasePDG.h>
#include <TParticle.h>
#include <TClonesArray.h>
#include <TClass.h>

ReadPar *sRPInstance;
STR      sRPFileName;
int      sNumEvents;
int      sRunType; 
int      sTables;
int      sModel;
int      sIntegrateSample;

ClassImp(TTherminator)

TTherminator::TTherminator():
  fCalka(0),
  fEvent(0),
  fPartDB(0)
{
  // Default constructor
  fPartDB = new ParticleDB();
}
TTherminator::TTherminator(const TTherminator & therm) :
  TGenerator(therm), 
  fCalka(0),
  fEvent(0),
  fPartDB(0)
{
  // Copy constructor
  //  fPartDB = new ParticleDB();
  if (fCalka) delete fCalka;
  fCalka = new Integrator(*therm.fCalka);
  if (fEvent) delete fEvent;
  fEvent = new Event(*therm.fEvent);
  if (fPartDB) delete fPartDB;
  fPartDB = new ParticleDB();
}
TTherminator& TTherminator::operator=(const TTherminator & therm)
{
  if (this != &therm) {
    fCalka = therm.fCalka;
    fEvent = therm.fEvent;
    delete fPartDB;
    fPartDB = new ParticleDB();
  }

  return *this;
}

TTherminator::~TTherminator()
{
  // Destructor
  if (sRPInstance) delete sRPInstance;
  if (fEvent) delete fEvent;
  if (fCalka) delete fCalka;
  if (fPartDB) delete fPartDB;
}

void TTherminator::ReadParameters()
{
  // Read in global model parameters
  STR tModel;
  STR tTable;
  
  try {
    sNumEvents = atoi(sRPInstance->getPar("NumberOfEvents").Data());
    tModel = sRPInstance->getPar("FreezeOutModel");
    if (tModel.Contains("SingleFreezeOut"))
      sModel = 0;
/*MCH begin*/
    else if (tModel.Contains("Lhyquid3D"))
      sModel = 10;
    else if (tModel.Contains("Lhyquid2D"))
      sModel = 11;
/*MCH end*/      
    else if (tModel.Contains("BlastWaveVTDelay"))
      sModel = 6;
    else if (tModel.Contains("BlastWaveVT"))
      sModel = 2;
    else if (tModel.Contains("BlastWaveVLinearFormation"))
      sModel = 7;
    else if (tModel.Contains("BlastWaveVLinearDelay"))
      sModel = 8;
    else if (tModel.Contains("BlastWaveVLinear"))
      sModel = 4;
    else if (tModel.Contains("FiniteHydro"))
      sModel = 5;
    else {
      PRINT_MESSAGE("Unknown model type: " << tModel.Data());
      PRINT_MESSAGE("Please provide the proper name");
      exit(0);
    }
    if (atoi(sRPInstance->getPar("Randomize").Data()) == 1)
      sRunType = 4;
    else
      sRunType = 3;
    tTable = sRPInstance->getPar("TableType");
    if (tTable.Contains("Mathematica"))
      sTables = 0;
    else if (tTable.Contains("SHARE"))
      sTables = 1;
    else {
      PRINT_MESSAGE("Unknown table type: " << tTable.Data());
      exit(0);
    }
    sIntegrateSample = atoi(sRPInstance->getPar("NumberOfIntegrateSamples").Data());
    fInputDir = sRPInstance->getPar("InputDirSHARE");
  }
  catch (STR tError) {
    PRINT_DEBUG_1("RunBFPW::ReadParameters - Caught exception " << tError);
    PRINT_MESSAGE("Did not find one of the neccessary parameters in the parameters file.");
    exit(0);
  }
}

void        TTherminator::Initialize(){
  // Initialize global variables
  // and particle and decay tables
  Parser       *tParser;
  
  sRPFileName = "therminator.in";
  sRPInstance = new ReadPar(sRPFileName.Data());

  ReadParameters();
  
  tParser = new Parser(fPartDB);

  tParser->ReadShare();
    
  fCalka = new Integrator(sIntegrateSample);
  fCalka->ReadMultInteg(fPartDB);


  // Read in particle table from share
  // and add missing particle type to TDatabasePDG

  char str[50];
  char str1[200];
    
  //  ParticleType *tPartBuf;

  TDatabasePDG *tInstance = TDatabasePDG::Instance();
  TParticlePDG *tParticleType;

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

  ifstream in((fInputDir+"/"+"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);
	    //	    printf("Added particle %s with PDG PID %d charge %d", str, (int) mc, charge);
	    //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();
    }

  delete tParser;
}

void        TTherminator::GenerateEvent()
{
  // Run single event generation
  if ((sRunType == 3) || (sRunType == 4))
    {
      if (sRunType==4)
	fCalka->Randomize();

      if (!fEvent) fEvent = new Event(fPartDB, fCalka);
      if (sRunType == 4) fEvent->Randomize();
      
      fEvent->GenerateEvent();
      fEvent->DecayParticles();
      //      fEvent->WriteEvent(0);
    }  
}

Int_t       TTherminator::ImportParticles(TClonesArray *particles, Option_t */*option*/)
{
  // Import particles from a generated event into an external array
  const double kFmToGev = 0.197327;

  if (particles == 0) return 0;
  TClonesArray &particlesR = *particles;
  particlesR.Clear();
  Int_t nump = 0;
  if (!fEvent) return 0;
  Int_t numpart = fEvent->GetParticleCount();
  printf("\n TTherminator: Therminator stack contains %d particles.\n", numpart);
  for (Int_t iPart=0; iPart<numpart; iPart++) {
    Particle *tPart = fEvent->GetParticleOfCount(iPart);
    Int_t tFather;
    tFather = tPart->GetFather();

    if (tFather> -1) {
      TParticle *mother = (TParticle*) (particlesR.UncheckedAt(tFather+1));	   
      mother->SetLastDaughter(iPart);
      if (mother->GetFirstDaughter()==-1)
 	mother->SetFirstDaughter(iPart);
    }
    nump++;
    //    printf("Putting %d %d %lf %d\n", tPart->GetParticleType()->GetPDGCode(), iPart, tPart->px, tPart);
    new (particlesR[iPart]) TParticle(tPart->GetParticleType()->GetPDGCode(), tPart->HadDecayed(),
 				      tFather, -1, -1, -1,
 				      tPart->px, tPart->py, tPart->pz, tPart->GetEnergy() ,
 				      tPart->rx*1.e-13*kFmToGev, tPart->ry*1.e-13*kFmToGev, tPart->rz*1.e-13*kFmToGev, tPart->rt*1.e-13*kFmToGev/3e10);
    particlesR[iPart]->SetUniqueID(iPart);
  }
  
  return nump;
}

TObjArray*  TTherminator::ImportParticles(Option_t */*option*/)
{
  // Import generated particles into an internal array
  const double kFmToGev = 0.197327;

  Int_t nump = 0;
  fParticles->Clear();
  if (!fEvent) return 0;
  Int_t numpart = fEvent->GetParticleCount();
  printf("\n TTherminator: Therminator stack contains %d particles.", numpart);
  fEvent->InitParticleScan();
  for (Int_t iPart=0; iPart<numpart; iPart++) {
    Particle *tPart = fEvent->GetNextParticle();
    Int_t tFather;
    tFather = tPart->GetFather();
    
    if (tFather> -1) {
      TParticle *mother = (TParticle*) (fParticles->UncheckedAt(tFather));	   
      mother->SetLastDaughter(iPart);
      if (mother->GetFirstDaughter()==-1)
 	mother->SetFirstDaughter(iPart);
    }
    TParticle* p = new TParticle(tPart->GetParticleType()->GetPDGCode(), tPart->HadDecayed(),
				 tFather, -1, -1, -1,
				 tPart->px, tPart->py, tPart->pz, tPart->GetEnergy() ,
				 tPart->rx*1.e-13*kFmToGev, tPart->ry*1.e-13*kFmToGev, tPart->rz*1.e-13*kFmToGev, tPart->rt*1.e-13*kFmToGev/3e10);
    p->SetUniqueID(iPart);
    fParticles->Add(p);
  }
  nump = fEvent->GetParticleCount();
  
  return fParticles;
  
}
Commit	Line	Data
2e967919	1	///////////////////////////////////////////////////////////////////////////////
	2	// //
	3	// TTherminator: global interface class to the Therminator model. //
	4	// Initialized global variables and runs the event generation //
	5	// //
	6	///////////////////////////////////////////////////////////////////////////////
	7	/******************************************************************************
	8	* T H E R M I N A T O R *
	9	* THERMal heavy-IoN generATOR *
	10	* version 1.0 *
	11	* *
	12	* Authors of the model: Wojciech Broniowski, Wojciech.Broniowski@ifj.edu.pl, *
	13	* Wojciech Florkowski, Wojciech.Florkowski@ifj.edu.pl *
	14	* Authors of the code: Adam Kisiel, kisiel@if.pw.edu.pl *
	15	* Tomasz Taluc, ttaluc@if.pw.edu.pl *
	16	* Code designers: Adam Kisiel, Tomasz Taluc, Wojciech Broniowski, *
	17	* Wojciech Florkowski *
	18	* *
	19	* For the detailed description of the program and furhter references *
	20	* to the description of the model plesase refer to: nucl-th/0504047, *
	21	* accessibile at: http://www.arxiv.org/nucl-th/0504047 *
	22	* *
	23	* Homepage: http://hirg.if.pw.edu.pl/en/therminator/ *
	24	* *
	25	* This code can be freely used and redistributed. However if you decide to *
	26	* make modifications to the code, please contact the authors, especially *
	27	* if you plan to publish the results obtained with such modified code. *
	28	* Any publication of results obtained using this code must include the *
	29	* reference to nucl-th/0504047 and the published version of it, when *
	30	* available. *
	31	* *
	32	*****************************************************************************/
	33	#include <TTherminator.h>
	34	#include <fstream>
	35	#include <TDatabasePDG.h>
	36	#include <TParticle.h>
	37	#include <TClonesArray.h>
0ff6d0b8	38	#include <TClass.h>
2e967919	39
	40	ReadPar *sRPInstance;
	41	STR sRPFileName;
	42	int sNumEvents;
	43	int sRunType;
	44	int sTables;
	45	int sModel;
	46	int sIntegrateSample;
	47
	48	ClassImp(TTherminator)
	49
	50	TTherminator::TTherminator():
	51	fCalka(0),
	52	fEvent(0),
	53	fPartDB(0)
	54	{
	55	// Default constructor
	56	fPartDB = new ParticleDB();
	57	}
	58	TTherminator::TTherminator(const TTherminator & therm) :
1aed17af	59	TGenerator(therm),
	60	fCalka(0),
	61	fEvent(0),
	62	fPartDB(0)
2e967919	63	{
2e967919	64	// Copy constructor
5c8cdc24	65	// fPartDB = new ParticleDB();
1aed17af	66	if (fCalka) delete fCalka;
	67	fCalka = new Integrator(*therm.fCalka);
	68	if (fEvent) delete fEvent;
	69	fEvent = new Event(*therm.fEvent);
	70	if (fPartDB) delete fPartDB;
	71	fPartDB = new ParticleDB();
5c8cdc24	72	}
	73	TTherminator& TTherminator::operator=(const TTherminator & therm)
	74	{
	75	if (this != &therm) {
	76	fCalka = therm.fCalka;
	77	fEvent = therm.fEvent;
	78	delete fPartDB;
	79	fPartDB = new ParticleDB();
	80	}
	81
	82	return *this;
2e967919	83	}
5c8cdc24	84
2e967919	85	TTherminator::~TTherminator()
	86	{
	87	// Destructor
	88	if (sRPInstance) delete sRPInstance;
	89	if (fEvent) delete fEvent;
	90	if (fCalka) delete fCalka;
	91	if (fPartDB) delete fPartDB;
	92	}
	93
	94	void TTherminator::ReadParameters()
	95	{
	96	// Read in global model parameters
	97	STR tModel;
	98	STR tTable;
	99
	100	try {
	101	sNumEvents = atoi(sRPInstance->getPar("NumberOfEvents").Data());
	102	tModel = sRPInstance->getPar("FreezeOutModel");
	103	if (tModel.Contains("SingleFreezeOut"))
	104	sModel = 0;
	105	/MCH begin/
	106	else if (tModel.Contains("Lhyquid3D"))
	107	sModel = 10;
	108	else if (tModel.Contains("Lhyquid2D"))
	109	sModel = 11;
	110	/MCH end/
	111	else if (tModel.Contains("BlastWaveVTDelay"))
	112	sModel = 6;
	113	else if (tModel.Contains("BlastWaveVT"))
	114	sModel = 2;
	115	else if (tModel.Contains("BlastWaveVLinearFormation"))
	116	sModel = 7;
	117	else if (tModel.Contains("BlastWaveVLinearDelay"))
	118	sModel = 8;
	119	else if (tModel.Contains("BlastWaveVLinear"))
	120	sModel = 4;
	121	else if (tModel.Contains("FiniteHydro"))
	122	sModel = 5;
	123	else {
	124	PRINT_MESSAGE("Unknown model type: " << tModel.Data());
	125	PRINT_MESSAGE("Please provide the proper name");
	126	exit(0);
	127	}
	128	if (atoi(sRPInstance->getPar("Randomize").Data()) == 1)
	129	sRunType = 4;
	130	else
	131	sRunType = 3;
	132	tTable = sRPInstance->getPar("TableType");
	133	if (tTable.Contains("Mathematica"))
	134	sTables = 0;
	135	else if (tTable.Contains("SHARE"))
	136	sTables = 1;
	137	else {
	138	PRINT_MESSAGE("Unknown table type: " << tTable.Data());
	139	exit(0);
	140	}
	141	sIntegrateSample = atoi(sRPInstance->getPar("NumberOfIntegrateSamples").Data());
ae89e57a	142	fInputDir = sRPInstance->getPar("InputDirSHARE");
2e967919	143	}
	144	catch (STR tError) {
	145	PRINT_DEBUG_1("RunBFPW::ReadParameters - Caught exception " << tError);
	146	PRINT_MESSAGE("Did not find one of the neccessary parameters in the parameters file.");
	147	exit(0);
	148	}
	149	}
	150
	151	void TTherminator::Initialize(){
	152	// Initialize global variables
	153	// and particle and decay tables
	154	Parser *tParser;
	155
	156	sRPFileName = "therminator.in";
	157	sRPInstance = new ReadPar(sRPFileName.Data());
	158
	159	ReadParameters();
	160
	161	tParser = new Parser(fPartDB);
	162
	163	tParser->ReadShare();
	164
	165	fCalka = new Integrator(sIntegrateSample);
	166	fCalka->ReadMultInteg(fPartDB);
	167
	168
	169	// Read in particle table from share
	170	// and add missing particle type to TDatabasePDG
	171
	172	char str[50];
	173	char str1[200];
	174
	175	// ParticleType *tPartBuf;
	176
	177	TDatabasePDG *tInstance = TDatabasePDG::Instance();
	178	TParticlePDG *tParticleType;
	179
	180	// AliWarning(Form("Reading particle types from particles.data"));
ae89e57a	181
	182	ifstream in((fInputDir+"/"+"particles.data").Data());
	183	// ifstream in("particles.data");
2e967919	184
	185	int charge;
	186
	187	int number=0;
	188	if ((in) && (in.is_open()))
	189	{
	190	//START OF HEAD-LINE
	191	in.ignore(200,'\n');
	192	in.ignore(200,'\n');
	193	in.ignore(200,'\n');
	194	//END OF HEAD-LINE
	195
	196	while (in>>str)
	197	{
	198	if (/(str == '#')\|\|/(str<65)\|\|(*str>122))
	199	{
	200	in.getline(str1,200);
	201	continue;
	202	}
	203	double mass, gamma, spin, tI3, tI, q, s, aq, as, c, ac, mc;
	204
	205	in>>mass>>gamma>>spin>>tI>>tI3>>q>>s>>aq>>as>>c>>ac>>mc;
	206	number++;
	207	tParticleType = tInstance->GetParticle((int) mc);
	208	if (!tParticleType) {
	209	charge = 0;
	210	if (strstr(str, "plu")) charge = 1;
	211	if (strstr(str, "min")) charge = -1;
	212	if (strstr(str, "plb")) charge = -1;
	213	if (strstr(str, "mnb")) charge = 1;
	214	if (strstr(str, "plp")) charge = 2;
	215	if (strstr(str, "ppb")) charge = -2;
	216	tInstance->AddParticle(str, str, mass, gamma == 0.0 ? 1:0, gamma, charge , "meson", (int) mc);
	217	// printf("Added particle %s with PDG PID %d charge %d", str, (int) mc, charge);
	218	//AliWarning(Form("Added particle %s with PDG PID %d charge %d", str, (int) mc, charge));
	219	// 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));
	220
	221	}
	222	}
	223	in.close();
	224	}
09b7c66c	225
09b7c66c	226	delete tParser;
2e967919	227	}
	228
	229	void TTherminator::GenerateEvent()
	230	{
	231	// Run single event generation
	232	if ((sRunType == 3) \|\| (sRunType == 4))
	233	{
	234	if (sRunType==4)
	235	fCalka->Randomize();
	236
	237	if (!fEvent) fEvent = new Event(fPartDB, fCalka);
	238	if (sRunType == 4) fEvent->Randomize();
	239
	240	fEvent->GenerateEvent();
	241	fEvent->DecayParticles();
	242	// fEvent->WriteEvent(0);
	243	}
	244	}
	245
c61a7285	246	Int_t TTherminator::ImportParticles(TClonesArray particles, Option_t /option/)
2e967919	247	{
2e967919	248	// Import particles from a generated event into an external array
c9c8970c	249	const double kFmToGev = 0.197327;
c9c8970c	250
2e967919	251	if (particles == 0) return 0;
	252	TClonesArray &particlesR = *particles;
	253	particlesR.Clear();
	254	Int_t nump = 0;
	255	if (!fEvent) return 0;
	256	Int_t numpart = fEvent->GetParticleCount();
0ff6d0b8	257	printf("\n TTherminator: Therminator stack contains %d particles.\n", numpart);
2e967919	258	for (Int_t iPart=0; iPart<numpart; iPart++) {
	259	Particle *tPart = fEvent->GetParticleOfCount(iPart);
	260	Int_t tFather;
	261	tFather = tPart->GetFather();
	262
	263	if (tFather> -1) {
	264	TParticle mother = (TParticle) (particlesR.UncheckedAt(tFather+1));
	265	mother->SetLastDaughter(iPart);
	266	if (mother->GetFirstDaughter()==-1)
	267	mother->SetFirstDaughter(iPart);
	268	}
	269	nump++;
	270	// printf("Putting %d %d %lf %d\n", tPart->GetParticleType()->GetPDGCode(), iPart, tPart->px, tPart);
	271	new (particlesR[iPart]) TParticle(tPart->GetParticleType()->GetPDGCode(), tPart->HadDecayed(),
0ff6d0b8	272	tFather, -1, -1, -1,
	273	tPart->px, tPart->py, tPart->pz, tPart->GetEnergy() ,
	274	tPart->rx1.e-13kFmToGev, tPart->ry1.e-13kFmToGev, tPart->rz1.e-13kFmToGev, tPart->rt1.e-13kFmToGev/3e10);
2e967919	275	particlesR[iPart]->SetUniqueID(iPart);
	276	}
	277
	278	return nump;
	279	}
	280
c61a7285	281	TObjArray* TTherminator::ImportParticles(Option_t /option*/)
2e967919	282	{
2e967919	283	// Import generated particles into an internal array
c9c8970c	284	const double kFmToGev = 0.197327;
c9c8970c	285
2e967919	286	Int_t nump = 0;
	287	fParticles->Clear();
	288	if (!fEvent) return 0;
	289	Int_t numpart = fEvent->GetParticleCount();
	290	printf("\n TTherminator: Therminator stack contains %d particles.", numpart);
	291	fEvent->InitParticleScan();
	292	for (Int_t iPart=0; iPart<numpart; iPart++) {
	293	Particle *tPart = fEvent->GetNextParticle();
	294	Int_t tFather;
	295	tFather = tPart->GetFather();
	296
	297	if (tFather> -1) {
	298	TParticle mother = (TParticle) (fParticles->UncheckedAt(tFather));
	299	mother->SetLastDaughter(iPart);
	300	if (mother->GetFirstDaughter()==-1)
	301	mother->SetFirstDaughter(iPart);
	302	}
	303	TParticle* p = new TParticle(tPart->GetParticleType()->GetPDGCode(), tPart->HadDecayed(),
	304	tFather, -1, -1, -1,
	305	tPart->px, tPart->py, tPart->pz, tPart->GetEnergy() ,
c9c8970c	306	tPart->rx1.e-13kFmToGev, tPart->ry1.e-13kFmToGev, tPart->rz1.e-13kFmToGev, tPart->rt1.e-13kFmToGev/3e10);
2e967919	307	p->SetUniqueID(iPart);
	308	fParticles->Add(p);
	309	}
	310	nump = fEvent->GetParticleCount();
	311
	312	return fParticles;
	313
	314	}