This commit was generated by cvs2svn to compensate for changes in r23244,
authorhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 8 Jan 2008 13:59:05 +0000 (13:59 +0000)
committerhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 8 Jan 2008 13:59:05 +0000 (13:59 +0000)
which included commits to RCS files with non-trunk default branches.

33 files changed:
TTherminator/AliGenTherminator.cxx [new file with mode: 0644]
TTherminator/AliGenTherminator.h [new file with mode: 0644]
TTherminator/Makefile [new file with mode: 0644]
TTherminator/PROOF-INF.TTherminator/BUILD.sh [new file with mode: 0755]
TTherminator/PROOF-INF.TTherminator/SETUP.C [new file with mode: 0644]
TTherminator/TTherminator.cxx [new file with mode: 0644]
TTherminator/TTherminator.h [new file with mode: 0644]
TTherminator/TTherminatorLinkDef.h [new file with mode: 0644]
TTherminator/Therminator/DecayChannel.cxx [new file with mode: 0644]
TTherminator/Therminator/DecayChannel.h [new file with mode: 0644]
TTherminator/Therminator/DecayTable.cxx [new file with mode: 0644]
TTherminator/Therminator/DecayTable.h [new file with mode: 0644]
TTherminator/Therminator/Event.cxx [new file with mode: 0644]
TTherminator/Therminator/Event.h [new file with mode: 0644]
TTherminator/Therminator/Hypersurface.cxx [new file with mode: 0644]
TTherminator/Therminator/Hypersurface.h [new file with mode: 0644]
TTherminator/Therminator/Integrator.cxx [new file with mode: 0644]
TTherminator/Therminator/Integrator.h [new file with mode: 0644]
TTherminator/Therminator/Makefile [new file with mode: 0644]
TTherminator/Therminator/Parser.cxx [new file with mode: 0644]
TTherminator/Therminator/Parser.h [new file with mode: 0644]
TTherminator/Therminator/Particle.cxx [new file with mode: 0644]
TTherminator/Therminator/Particle.h [new file with mode: 0644]
TTherminator/Therminator/ParticleDB.cxx [new file with mode: 0644]
TTherminator/Therminator/ParticleDB.h [new file with mode: 0644]
TTherminator/Therminator/ParticleDecayer.cxx [new file with mode: 0644]
TTherminator/Therminator/ParticleDecayer.h [new file with mode: 0644]
TTherminator/Therminator/ParticleType.cxx [new file with mode: 0644]
TTherminator/Therminator/ParticleType.h [new file with mode: 0644]
TTherminator/Therminator/ReadPar.cxx [new file with mode: 0644]
TTherminator/Therminator/ReadPar.h [new file with mode: 0644]
TTherminator/Therminator/THGlobal.h [new file with mode: 0644]
TTherminator/libTTherminator.pkg [new file with mode: 0644]

diff --git a/TTherminator/AliGenTherminator.cxx b/TTherminator/AliGenTherminator.cxx
new file mode 100644 (file)
index 0000000..dffe99c
--- /dev/null
@@ -0,0 +1,261 @@
+// 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;
+}
diff --git a/TTherminator/AliGenTherminator.h b/TTherminator/AliGenTherminator.h
new file mode 100644 (file)
index 0000000..48d025d
--- /dev/null
@@ -0,0 +1,97 @@
+#ifndef ALIGENTHERMINATOR_H
+#define ALIGENTHERMINATOR_H
+
+// ALICE event generator based on the THERMINATOR model
+// It invokes the event generation and then puts the resulting
+// event on the stack. All the parameters of the model are 
+// accessible via this interface
+// Author: Adam.Kisiel@cern.ch
+
+class TF1;
+class TClonesArray;
+
+#include "AliGenerator.h"
+#include "AliDecayer.h"
+#include "AliGenMC.h"
+#include "TTherminator.h"
+
+class AliGenTherminator : public AliGenMC
+{
+ public:
+
+  AliGenTherminator();
+  AliGenTherminator(Int_t npart);
+  virtual ~AliGenTherminator();
+  virtual void Generate();
+  virtual void Init();
+  // Set model parameters
+
+  void SetFileName(const char *infilename);
+  void SetEventNumberInFile(int evnum);
+
+  void SetTemperature(Double_t temp);       
+  void SetMiuI(Double_t miu);              
+  void SetMiuS(Double_t miu);              
+  void SetMiuB(Double_t miu);              
+  void SetAlfaRange(Double_t range);         
+  void SetRapRange(Double_t range);          
+
+  void SetRhoMax(Double_t rho);  
+  void SetTau(Double_t tau);     
+  void SetBWA(Double_t bwa);     
+  void SetBWVt(Double_t bwv);    
+  void SetBWDelay(Double_t bwd); 
+
+  void SetModel(const char *model);
+
+ protected:
+  void     ReadShareParticleTable();
+  void     CreateTherminatorInputFile();
+
+  Int_t    fNt;                // CurrentTrack;
+  Int_t    fEventNumber;       // Number of the event to read
+  TString  fFileName;          // FileName of the file with events
+
+  TString  fFreezeOutModel;    // FreezeOut model to use
+
+  // Parameters common for all models
+  Double_t fTemperature;       // Freeze-out temperature [GeV/c]
+  Double_t fMiuI;              // Isospin chemical potential
+  Double_t fMiuS;              // Strance chemical potential
+  Double_t fMiuB;              // Baryonic chemical potential
+  Double_t fAlfaRange;         // Spatial rapidity range for primordial particles
+  Double_t fRapRange;          // Momentum rapidity range for primordial particles
+  
+  // Model dependent paramters
+  Double_t fRhoMax;            // Transverse bound of the system
+  Double_t fTau;               // Proper freeze-out time parameter
+  Double_t fBWA;               // Blast-wave freeze-out surface slope in t-rho plane
+  Double_t fBWVt;              // Blast-wave velocity profile parameter
+  Double_t fBWDelay;           // Blast-wave time delay parameter
+
+ private:
+
+  ClassDef(AliGenTherminator,1) // Hijing parametrisation generator
+};
+#endif
+
+inline void AliGenTherminator::SetTemperature(Double_t temp) { fTemperature = temp; }
+inline void AliGenTherminator::SetMiuI(Double_t miu) { fMiuI = miu; }
+inline void AliGenTherminator::SetMiuS(Double_t miu) { fMiuS = miu; }
+inline void AliGenTherminator::SetMiuB(Double_t miu) { fMiuB = miu; }
+inline void AliGenTherminator::SetAlfaRange(Double_t range) { fAlfaRange = range; }
+inline void AliGenTherminator::SetRapRange(Double_t range) { fRapRange = range; }
+
+inline void AliGenTherminator::SetRhoMax(Double_t rho) { fRhoMax = rho; }
+inline void AliGenTherminator::SetTau(Double_t tau) { fTau = tau; }
+inline void AliGenTherminator::SetBWA(Double_t bwa) { fBWA = bwa; }
+inline void AliGenTherminator::SetBWVt(Double_t bwv) { fBWVt = bwv; }
+inline void AliGenTherminator::SetBWDelay(Double_t bwd) { fBWDelay = bwd; }
+
+
+
+
+
+
+
+
diff --git a/TTherminator/Makefile b/TTherminator/Makefile
new file mode 100644 (file)
index 0000000..f14f54e
--- /dev/null
@@ -0,0 +1,86 @@
+# $Id$
+
+include $(ROOTSYS)/test/Makefile.arch
+
+default-target: libTTherminator.so
+
+ALICEINC      = -I. -I$(ALICE_ROOT)/include -ITherminator
+
+### define include dir for local case and par case
+ifneq ($(ALICE_ROOT),)
+  ALICEINC += -I$(ALICE_ROOT)/include -I$(ALICE_ROOT)/PYTHIA6 -I$(ALICE_ROOT)/EVGEN -I$(ALICE_ROOT)/RAW -I$(ALICE_ROOT)/TPC
+else
+  ifneq ($(STEERBase_INCLUDE),)
+    ALICEINC += -I../$(STEERBase_INCLUDE)
+  endif
+  ifneq ($(ESD_INCLUDE),)
+    ALICEINC += -I../$(ESD_INCLUDE)
+  endif
+endif
+
+PACKAGE = TTherminator
+include lib$(PACKAGE).pkg
+
+DHDR_TTherminator := $(DHDR)
+HDRS_TTherminator := $(HDRS)
+SRCS_TTherminator := $(SRCS) G__$(PACKAGE).cxx
+OBJS_TTherminator := $(SRCS_TTherminator:.cxx=.o)
+
+PARFILE       = $(PACKAGE).par
+
+lib$(PACKAGE).so: $(OBJS_TTherminator)
+       @echo "Linking" $@ ...
+       @/bin/rm -f $@
+ifeq ($(PLATFORM),macosx)
+       @$(LD) -bundle -undefined $(UNDEFOPT) $(LDFLAGS) $^ -o $@
+else
+       @$(LD) $(SOFLAGS) $(LDFLAGS) $^ -o $@
+endif
+       @chmod a+x $@
+       @echo "done"
+
+%.o:    %.cxx %.h
+       $(CXX) $(CXXFLAGS) -c $< -o $@ $(ALICEINC)
+
+G__TTherminator.cxx G__TTherminator.h: $(HDRS_TTherminator) $(DHDR_TTherminator)
+       @echo "Generating dictionary ..."
+       rootcint -f $@ -c $(ALICEINC) $^
+
+clean:
+       @rm -f $(OBJS_TTherminator) *.so G__TTherminator.* $(PARFILE)
+
+### CREATE PAR FILE
+
+$(PARFILE): $(patsubst %,$(PACKAGE)/%,$(filter-out G__%, $(HDRS_TTherminator) $(SRCS_TTherminator) $(DHDR_TTherminator) Makefile Makefile.arch lib$(PACKAGE).pkg PROOF-INF))
+       @echo "Creating archive" $@ ...
+       @tar cfzh $@ $(PACKAGE)
+       @rm -rf $(PACKAGE)
+       @echo "done"
+
+$(PACKAGE)/Makefile: Makefile #.$(PACKAGE)
+       @echo Copying $< to $@ with transformations
+       @[ -d $(dir $@) ] || mkdir -p $(dir $@)
+       @sed 's/include \$$(ROOTSYS)\/test\/Makefile.arch/include Makefile.arch/' < $^ > $@
+
+$(PACKAGE)/Makefile.arch: $(ROOTSYS)/test/Makefile.arch
+       @echo Copying $< to $@
+       @[ -d $(dir $@) ] || mkdir -p $(dir $@)
+       @cp -a $^ $@
+
+$(PACKAGE)/PROOF-INF: PROOF-INF.$(PACKAGE)
+       @echo Copying $< to $@
+       @[ -d $(dir $@) ] || mkdir -p $(dir $@)
+       @cp -a -r $^ $@
+
+$(PACKAGE)/%: %
+       @echo Copying $< to $@
+       @[ -d $(dir $@) ] || mkdir -p $(dir $@)
+       @cp -a $< $@
+
+test-%.par: %.par
+       @echo "INFO: The file $< is now tested, in case of an error check in par-tmp."
+       @mkdir -p par-tmp
+       @cd par-tmp; tar xfz ../$<;     cd $(subst .par,,$<); PROOF-INF/BUILD.sh
+       @rm -rf par-tmp
+       @echo "INFO: Testing succeeded (already cleaned up)"
+
diff --git a/TTherminator/PROOF-INF.TTherminator/BUILD.sh b/TTherminator/PROOF-INF.TTherminator/BUILD.sh
new file mode 100755 (executable)
index 0000000..92b2328
--- /dev/null
@@ -0,0 +1,5 @@
+#! /bin/sh
+
+touch libTTherminator.pkg
+
+make libTTherminator.so
diff --git a/TTherminator/PROOF-INF.TTherminator/SETUP.C b/TTherminator/PROOF-INF.TTherminator/SETUP.C
new file mode 100644 (file)
index 0000000..27a5957
--- /dev/null
@@ -0,0 +1,18 @@
+void SETUP() {
+  CheckLoadLibrary("libTTherminator");
+
+  // Set the include paths
+  gROOT->ProcessLine(".include TTherminator");
+
+  // Set our location, so that other packages can find us
+  gSystem->Setenv("TTherminator_INCLUDE", "TTherminator");
+}
+
+Int_t CheckLoadLibrary(const char* library) {
+  // checks if a library is already loaded, if not loads the library
+
+  if (strlen(gSystem->GetLibraries(Form("%s.so", library), "", kFALSE)) > 0)
+    return 1;
+  
+  return gSystem->Load(library);
+}
diff --git a/TTherminator/TTherminator.cxx b/TTherminator/TTherminator.cxx
new file mode 100644 (file)
index 0000000..e256eb9
--- /dev/null
@@ -0,0 +1,287 @@
+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+// 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>
+
+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) :
+  fCalka(0),
+  fEvent(0),
+  fPartDB(0)
+{
+  // Copy constructor
+  fPartDB = new ParticleDB();
+}
+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());
+  }
+  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("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();
+    }
+}
+
+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
+  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.", 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, tPart->ry*1.e-13, tPart->rz*1.e-13, tPart->rt*1.e-13/300000000 
+                                     );
+    particlesR[iPart]->SetUniqueID(iPart);
+  }
+  
+  return nump;
+}
+
+TObjArray*  TTherminator::ImportParticles(Option_t *option)
+{
+  // Import generated particles into an internal array
+  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, tPart->ry*1.e-13, tPart->rz*1.e-13, tPart->rt*1.e-13/300000000 
+                                );
+    p->SetUniqueID(iPart);
+    fParticles->Add(p);
+  }
+  nump = fEvent->GetParticleCount();
+  
+  return fParticles;
+  
+}
diff --git a/TTherminator/TTherminator.h b/TTherminator/TTherminator.h
new file mode 100644 (file)
index 0000000..881564f
--- /dev/null
@@ -0,0 +1,68 @@
+///////////////////////////////////////////////////////////////////////////////
+// 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.                                                                 *
+ *                                                                            *
+ *****************************************************************************/
+#ifndef TTHERMINATOR_H
+#define TTHERMINATOR_H
+
+#include <TGenerator.h>
+#include <map>
+#include <THGlobal.h>
+#include "ReadPar.h"
+#include "Parser.h"
+#include "ParticleDB.h"
+#include "ParticleType.h"
+#include "DecayTable.h"
+#include "Event.h"
+
+class TTherminator: public TGenerator {
+ public:
+  TTherminator();
+  TTherminator(const TTherminator & therm);
+  virtual ~TTherminator();
+  
+  virtual void        ReadParameters();
+  
+  virtual void        Initialize();
+  
+  virtual void        GenerateEvent();
+  
+  virtual Int_t       ImportParticles(TClonesArray *particles, Option_t *option="");
+  virtual TObjArray*  ImportParticles(Option_t *option="");
+
+ private:
+  int Run();
+  Integrator *fCalka; // Integrator class
+  Event      *fEvent; // The therminator event
+  ParticleDB *fPartDB;// Particle properties database
+
+  ClassDef(TTherminator,1) // Hijing parametrisation generator
+};
+
+#endif
diff --git a/TTherminator/TTherminatorLinkDef.h b/TTherminator/TTherminatorLinkDef.h
new file mode 100644 (file)
index 0000000..7768e23
--- /dev/null
@@ -0,0 +1,9 @@
+#ifdef __CINT__
+#pragma link off all globals;
+#pragma link off all classes;
+#pragma link off all functions;
+#pragma link C++ class TTherminator+;
+#pragma link C++ class AliGenTherminator+;
+#endif
diff --git a/TTherminator/Therminator/DecayChannel.cxx b/TTherminator/Therminator/DecayChannel.cxx
new file mode 100644 (file)
index 0000000..de8becd
--- /dev/null
@@ -0,0 +1,101 @@
+/******************************************************************************
+ *                      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 "DecayChannel.h"
+
+DecayChannel::DecayChannel() :
+  mParticleType1(0), mParticleType2(0), mParticleType3(-1), mBranchRatio(0.0)
+{
+}
+
+DecayChannel::DecayChannel(const DecayChannel& aChannel) 
+{
+  mBranchRatio = aChannel.GetBranchingRatio();
+  mParticleType1 = aChannel.GetParticle1();
+  mParticleType2 = aChannel.GetParticle2();
+  mParticleType3 = aChannel.GetParticle3();
+}
+
+DecayChannel::DecayChannel(double aBranchRatio, int aPartType1, int aPartType2, int aPartType3) :
+  mParticleType1(aPartType1), mParticleType2(aPartType2), mParticleType3(aPartType3), mBranchRatio(aBranchRatio)
+{
+}
+
+DecayChannel::~DecayChannel()
+{
+}
+
+int    
+DecayChannel::GetParticle1() const
+{
+  return mParticleType1;
+}
+
+int    
+DecayChannel::GetParticle2() const
+{
+  return mParticleType2;
+}
+
+int    
+DecayChannel::GetParticle3() const
+{
+  return mParticleType3;
+}
+
+double 
+DecayChannel::GetBranchingRatio() const
+{
+  return mBranchRatio;
+}
+
+int    
+DecayChannel::Is3Particle() const
+{
+  return (mParticleType3 != -1);
+}
+
+
+void   DecayChannel::SetParticle1(int aPartType1)
+{
+  mParticleType1 = aPartType1;
+}
+
+void   DecayChannel::SetParticle2(int aPartType2)
+{
+  mParticleType2 = aPartType2;
+}
+
+void   
+DecayChannel::SetParticle3(int aPartType3)
+{
+  mParticleType3 = aPartType3;
+}
+
+void   DecayChannel::SetBranchingRatio(double aRatio)
+{
+  mBranchRatio = aRatio;
+}
diff --git a/TTherminator/Therminator/DecayChannel.h b/TTherminator/Therminator/DecayChannel.h
new file mode 100644 (file)
index 0000000..0ca744e
--- /dev/null
@@ -0,0 +1,59 @@
+/******************************************************************************
+ *                      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.                                                                 *
+ *                                                                            *
+ *****************************************************************************/
+
+#ifndef _BFPW_DECAY_CHANNEL_
+#define _BFPW_DECAY_CHANNEL_
+#include "ParticleType.h"
+
+class DecayChannel
+{
+ public:
+  DecayChannel();
+  DecayChannel(double aBranchRatio, int aPartType1, int aPartType2, int aPartType3 = -1);
+  DecayChannel(const DecayChannel& aChannel);
+  ~DecayChannel();
+  
+  int    GetParticle1() const;
+  int    GetParticle2() const;
+  int    GetParticle3() const;
+  double GetBranchingRatio() const;
+  int    Is3Particle() const;
+  
+  void   SetParticle1(int aPartType1);
+  void   SetParticle2(int aPartType2);
+  void   SetParticle3(int aPartType3);
+  void   SetBranchingRatio(double aRatio);
+  
+ protected:
+  int    mParticleType1;
+  int    mParticleType2;
+  int    mParticleType3;
+  double mBranchRatio;
+  
+};
+
+#endif
diff --git a/TTherminator/Therminator/DecayTable.cxx b/TTherminator/Therminator/DecayTable.cxx
new file mode 100644 (file)
index 0000000..53bb773
--- /dev/null
@@ -0,0 +1,118 @@
+/******************************************************************************
+ *                      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 "DecayTable.h"
+
+using namespace std;
+
+DecayTable::DecayTable()
+{
+  mDecayChannels.clear();
+  mBranchingRatios.clear();
+}
+
+DecayTable::DecayTable(const DecayTable& aTable)
+{
+  
+  mDecayChannels.clear();
+  mBranchingRatios.clear();
+  for (int iter=0; iter<aTable.GetChannelCount(); iter++)
+    AddDecayChannel(*(aTable.GetDecayChannel(iter)));
+}
+
+DecayTable::~DecayTable()
+{
+}
+
+void 
+DecayTable::AddDecayChannel(DecayChannel aChannel)
+{
+  mDecayChannels.push_back(aChannel);
+  RecalculateBranchingRatios();
+}
+
+void 
+DecayTable::RecalculateBranchingRatios()
+{
+  float tSumRatio = 0.0;
+  float tCurRatio = 0.0;
+  
+  for (unsigned int tIter=0; tIter<mDecayChannels.size(); tIter++)
+    tSumRatio += mDecayChannels[tIter].GetBranchingRatio();
+  
+  for (unsigned int tIter=0; tIter<mDecayChannels.size(); tIter++)
+    {
+      tCurRatio += mDecayChannels[tIter].GetBranchingRatio()/tSumRatio;
+      if (mBranchingRatios.size() <= tIter)
+       mBranchingRatios.push_back(tCurRatio);
+      else
+       mBranchingRatios[tIter] = tCurRatio;
+    }
+}
+
+int  
+DecayTable::GetChannelCount() const
+{
+  return mDecayChannels.size()-1;
+}
+
+
+const DecayChannel* 
+DecayTable::GetDecayChannel(int aIndex) const
+{
+  return &(mDecayChannels[aIndex]);
+}
+
+float 
+DecayTable::GetDecayStep(int aIndex)
+{
+  return mBranchingRatios[aIndex];
+}
+
+int   
+DecayTable::ChooseDecayChannel(double aProb)
+{
+  unsigned int tChanIndex = 0;
+  while ((mBranchingRatios[tChanIndex] < aProb) && (tChanIndex < mDecayChannels.size()))
+    tChanIndex++;
+  
+  return tChanIndex;
+}
+
+int   
+DecayTable::ChooseDecayChannelOrNot(double aProb)
+{
+  float tSumRatio = 0.0;
+  
+  for (unsigned int tIter=0; tIter<mDecayChannels.size(); tIter++) {
+    if ((aProb > tSumRatio) && (aProb <= tSumRatio+mDecayChannels[tIter].GetBranchingRatio()))
+      return tIter;
+    tSumRatio += mDecayChannels[tIter].GetBranchingRatio();
+  }
+
+  return -1;
+}
+
diff --git a/TTherminator/Therminator/DecayTable.h b/TTherminator/Therminator/DecayTable.h
new file mode 100644 (file)
index 0000000..3272631
--- /dev/null
@@ -0,0 +1,58 @@
+/******************************************************************************
+ *                      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.                                                                 *
+ *                                                                            *
+ *****************************************************************************/
+#ifndef _BFPW_DECAY_TABLE_
+#define _BFPW_DECAY_TABLE_
+
+#include <vector>
+#include "THGlobal.h"
+#include "DecayChannel.h"
+
+class DecayTable
+{
+ public:
+  DecayTable();
+  DecayTable(const DecayTable& aTable);
+  ~DecayTable();
+  
+  void AddDecayChannel(DecayChannel aChannel);
+
+  int  GetChannelCount() const;
+  const DecayChannel* GetDecayChannel(int aIndex) const;
+  float GetDecayStep(int aIndex);
+  int   ChooseDecayChannel(double aProb);
+  int   ChooseDecayChannelOrNot(double aProb);
+  
+ private:
+  std::vector<DecayChannel> mDecayChannels;
+  std::vector<float>        mBranchingRatios;
+  
+  void RecalculateBranchingRatios();
+};
+
+
+
+#endif
diff --git a/TTherminator/Therminator/Event.cxx b/TTherminator/Therminator/Event.cxx
new file mode 100644 (file)
index 0000000..b816b93
--- /dev/null
@@ -0,0 +1,359 @@
+/******************************************************************************
+ *                      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 <iostream>
+#include <TFile.h>
+#include <TH1D.h>
+#include "Particle.h"
+#include "Event.h"
+#include "Integrator.h"
+#include "ParticleDecayer.h"
+#include "DecayTable.h"
+#include "ReadPar.h"
+
+extern ReadPar* sRPInstance;
+
+Event::Event(ParticleDB *aDB, Integrator* aInteg) :
+  mCurIter(NULL),
+  mScanStarted(0)
+{
+  mDB = aDB;
+  mInteg = aInteg;
+  mRandom = new TRandom2();
+#ifdef _ROOT_4_
+  mRandom->SetSeed2(31851, 14327);
+#else
+  mRandom->SetSeed(31851);
+#endif
+  mAverageMultiplicities.resize(mDB->GetParticleTypeCount());
+  mMultiplicities.resize(mDB->GetParticleTypeCount());
+  ReadParameters();
+  ReadMultiplicities();
+
+
+}
+
+Event::~Event()
+{
+  
+  mParticles.clear();
+  delete mRandom;
+}
+
+void 
+Event::GenerateEvent(int aSeed)
+{
+  Particle **tPartBuf;
+
+  mParticles.clear();
+#ifdef _ROOT_4_
+  if (aSeed) mRandom->SetSeed2(aSeed, (aSeed*2) % (7*11*23*31));
+#else
+  if (aSeed) mRandom->SetSeed(aSeed);
+#endif
+  GenerateMultiplicities();
+  for (unsigned int tIter=0; tIter<mAverageMultiplicities.size(); tIter++) {
+    mInteg->Generate(mDB->GetParticleType(tIter), mMultiplicities[tIter], &tPartBuf);
+    for (int tBufIter=0; tBufIter<mMultiplicities[tIter]; tBufIter++) {
+      AddParticle(tPartBuf[tBufIter]);
+      delete tPartBuf[tBufIter];
+    }
+    free(tPartBuf);
+  }
+}
+
+int  
+Event::GetParticleCount()
+{
+  return mParticles.size();
+}
+
+Particle* 
+Event::GetNextParticle()
+{
+  if (mScanStarted)
+    InitParticleScan();
+  if (mCurIter == mParticles.end())
+    return 0;
+
+  Particle *tBuf;
+  tBuf = &(*mCurIter);
+  mCurIter++;
+  
+  return tBuf;
+}
+
+void 
+Event::InitParticleScan()
+{
+  mCurIter = mParticles.begin();
+  mScanStarted = 1;
+}
+
+Particle* 
+Event::GetParticleOfCount(int aCount)
+{
+  if (aCount == 0) {
+    InitParticleScan();
+  }
+  else if (aCount == 1) {
+    InitParticleScan();
+    mCurIter++;
+  }
+  else {
+    mCurIter--;
+    mCurIter++;
+    mCurIter++;
+  }
+
+  if (mCurIter == mParticles.end())
+    return 0;
+  
+  Particle *tBuf;
+  tBuf = &(*mCurIter);
+  return tBuf;
+}
+
+void 
+Event::WriteEvent(int nEvent)
+{
+  ParticleListIter tPLIter;
+  tPLIter = mParticles.begin();
+  int tCount = 0;
+  ofstream os;
+  
+  if (nEvent == 0) {
+    try {
+      os.open(sRPInstance->getPar("EventOutputFile").Data());
+    }
+    catch (STR e) {
+      PRINT_MESSAGE("Very strange: No event output filename ???");
+      exit(0);
+    }
+  }
+  else {
+    try {
+      os.open(sRPInstance->getPar("EventOutputFile").Data(), ios::app);
+    }
+    catch (STR e) {
+      PRINT_MESSAGE("Very strange: No event output filename ???");
+      exit(0);
+    }
+  }
+    
+  
+
+  if (nEvent == 0) 
+    {
+      os << "Therminator text output\nall_particles_p_x\nTherminator " << endl;
+    }
+  os << nEvent+1 << "  " << GetParticleCount() << "    " << "0.000" << "       " << "0.000" << endl;
+
+  while (tPLIter != mParticles.end())
+    {
+      os.flags(ios::right);
+      os.precision(6);
+      os.width(6);
+      os << tCount << "   ";
+      (*tPLIter).WriteParticle(&os);
+      os << endl;
+
+      tCount++;
+      tPLIter++;
+    }
+}
+
+  
+void 
+Event::AddParticle(Particle* aParticle)
+{
+  if (0) {
+    ParticleListIter tPLIter;
+    tPLIter = mParticles.begin();
+    while (tPLIter != mParticles.end())
+      {
+       if (((*tPLIter).GetMass() == aParticle->GetMass()) &&
+           ((*tPLIter).GetI3() == aParticle->GetI3()) &&
+           ((*tPLIter).GetBarionN() == aParticle->GetBarionN()) &&
+           ((*tPLIter).GetStrangeness() == aParticle->GetStrangeness()))
+         break;
+       if ((*tPLIter).GetMass() < aParticle->GetMass())
+         break;
+       else if ((*tPLIter).GetMass() == aParticle->GetMass())
+         if ((*tPLIter).GetI3() < aParticle->GetI3())
+           break;
+         else if ((*tPLIter).GetI3() == aParticle->GetI3())
+           if ((*tPLIter).GetBarionN() < aParticle->GetBarionN())
+             break;
+           else if ((*tPLIter).GetBarionN() == aParticle->GetBarionN())
+             if ((*tPLIter).GetStrangeness() < aParticle->GetStrangeness())
+               break;
+       tPLIter++;
+      }
+    mParticles.insert(tPLIter, *aParticle);
+  }
+  else
+    mParticles.push_back(*aParticle);
+}
+
+void           
+Event::ReadMultiplicities()
+{
+  char   tName[200];
+  double tMult;
+  
+  char *tHash;
+  char  tMultName[100];
+  ifstream *fin = NULL;
+  
+  tHash = mInteg->ParameterHash();
+  
+  strcpy(tMultName, "fmultiplicity_");
+  strcat(tMultName, tHash);
+  strcat(tMultName, ".txt");
+
+  fin = new ifstream(tMultName);
+  if ((fin) && (fin->is_open())) {
+    PRINT_MESSAGE("Reading Multiplicities from " << tMultName);
+
+    while (!fin->eof())
+    {
+      (*fin) >> tName >> tMult;
+      PRINT_DEBUG_2(tName << " " <<  mDB->GetParticleTypeIndex(strdup(tName)) << " " << tMult);
+      mAverageMultiplicities[mDB->GetParticleTypeIndex(strdup(tName))] = tMult;
+    }
+    fin->close();
+  }
+}
+
+void           
+Event::GenerateMultiplicities()
+{
+  for (unsigned int tIter=0; tIter<mAverageMultiplicities.size(); tIter++)
+    mMultiplicities[tIter] = mRandom->Poisson(mAverageMultiplicities[tIter]);
+}
+
+void 
+Event::DecayParticles()
+{
+  Particle *tPart1, *tPart2, *tPart3, *tFather;
+  ParticleDecayer* tDecayer;
+  int tCount = 0;
+  
+  tDecayer = new ParticleDecayer(mDB);
+  tDecayer->SeedSet(mRandom->Integer(100000000));
+  
+  //  InitParticleScan();
+  while ((tFather = GetParticleOfCount(tCount))) 
+    {
+      if (tFather->GetParticleType()->GetGamma() >= 0.0)
+       if ((tFather->GetParticleType()->GetTable()) && (((DecayTable *) tFather->GetParticleType()->GetTable())->GetChannelCount()+1 > 0))
+         {
+           tDecayer->DecayParticle(tFather, &tPart1, &tPart2, &tPart3);
+#ifndef _RESCALE_CHANNELS_
+           if (!tPart1) {
+             tCount++;
+             continue;
+           }
+#endif
+#ifdef _NO_THREE_BODY_DECAYS_
+           if (tPart3){
+             tCount++;
+             continue;
+           }
+#endif
+#ifdef _OMIT_TWO_BODY_
+           if (tPart1 && tPart2 && !tPart3){
+             tCount++;
+             continue;
+           }
+#endif
+
+           tPart1->SetFather(tCount);
+           tPart2->SetFather(tCount);
+           if (tPart3)
+             tPart3->SetFather(tCount);
+           
+           AddParticle(tPart1);
+           AddParticle(tPart2);
+           
+           delete tPart1;
+           delete tPart2;
+
+           if (tPart3) {
+             AddParticle(tPart3);
+             delete tPart3;
+           }
+         }
+       else
+         {
+         }
+      tCount++;
+    }
+
+  delete tDecayer;
+}
+
+void 
+Event::Randomize()
+{
+  TDatime dat;
+  
+#ifdef _ROOT_4_
+  mRandom->SetSeed2(dat.Get() / 2 * 3, dat.Get() / 11 * 9);
+#else
+  mRandom->SetSeed(dat.Get() / 2 * 3);
+#endif
+}
+
+void           
+Event::ReadParameters()
+{
+  STR tEvFile;
+  STR tUseNegativeBinomial;
+
+  try {
+    tEvFile = sRPInstance->getPar("EventOutputFile");
+  }
+  catch (STR e) {
+    PRINT_DEBUG_1("Event::ReadParameters - Caught exception " << e);
+    PRINT_MESSAGE("Did not find one of the neccessary parameters in the parameters file.");
+    exit(0);
+  }
+  try {
+    tUseNegativeBinomial = sRPInstance->getPar("MultiplicityDistribution");
+    if (tUseNegativeBinomial.Contains("NegativeBinomial"))
+      mNegBin = 1;
+    else
+      mNegBin = 0;
+  }
+  catch (STR e) {
+    PRINT_MESSAGE("Using default multiplicty distribution: Poissonian");
+    mNegBin = 0;
+  }
+       
+}
+
diff --git a/TTherminator/Therminator/Event.h b/TTherminator/Therminator/Event.h
new file mode 100644 (file)
index 0000000..34d19cc
--- /dev/null
@@ -0,0 +1,78 @@
+/******************************************************************************
+ *                      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.                                                                 *
+ *                                                                            *
+ *****************************************************************************/
+#ifndef _BFPW_EVENT_
+#define _BFPW_EVENT_
+
+#include <vector>
+#include <list>
+#include <TRandom2.h>
+#include "THGlobal.h"
+#include "Particle.h"
+#include "ParticleDB.h"
+#include "Integrator.h"
+
+using namespace std;
+
+typedef list<Particle> ParticleList;
+typedef list<Particle>::iterator ParticleListIter;
+
+class Particle;
+
+class Event
+{
+ public:
+  Event(ParticleDB *aDB, Integrator *aInteg);
+  ~Event();
+
+  void GenerateEvent(int aSeed=0);
+  void DecayParticles();
+  
+  int  GetParticleCount();
+  void InitParticleScan();
+  Particle* GetNextParticle();
+  Particle* GetParticleOfCount(int aCount);
+  void WriteEvent(int nEvent=0);
+  void AddParticle(Particle* aParticle);
+  void Randomize();
+  
+ private:
+  ParticleList   mParticles;
+  vector<int>    mMultiplicities;
+  vector<double> mAverageMultiplicities;
+  void           GenerateMultiplicities();
+  void           ReadMultiplicities();
+  void           ReadParameters();
+  ParticleDB    *mDB;
+  Integrator    *mInteg;
+  TRandom2      *mRandom;
+  ParticleListIter mCurIter;
+  Int_t          mNegBin;
+  Int_t          mScanStarted;
+};
+
+
+#endif
diff --git a/TTherminator/Therminator/Hypersurface.cxx b/TTherminator/Therminator/Hypersurface.cxx
new file mode 100644 (file)
index 0000000..056df9c
--- /dev/null
@@ -0,0 +1,196 @@
+#include "Hypersurface.h"
+
+Hypersurface::Hypersurface(void) {
+// Read FOHSI.txt file and init parameters
+  FName = "fohsi.txt";
+  HSFile = new ifstream;
+  HSFile->open(FName);
+  if(HSFile->is_open()) {
+    HSFile->seekg (0, ios::beg);
+// phi
+    HSFile->getline(buff,100);
+    Np=atoi(buff);
+    HSFile->getline(buff,100);
+    ip=atof(buff);
+    HSFile->getline(buff,100);
+    fp=atof(buff);
+// zeta
+    HSFile->getline(buff,100);
+    Nz=atoi(buff);
+    HSFile->getline(buff,100);
+    iz=atof(buff);
+    HSFile->getline(buff,100);
+    fz=atof(buff);
+// freeze-out temperature
+    HSFile->getline(buff,100);
+    TFO=atof(buff);
+
+    HSFile->close();
+    dp = (fp-ip)/(Np-1);
+    dz = (fz-iz)/(Nz-1);
+
+// create 2D array
+    aArr   = new double* [Np];
+    vArr   = new double* [Np];
+    dArr   = new double* [Np];
+    DpdArr = new double* [Np];
+    DzdArr = new double* [Np];
+    for(i=0;i<Np;i++) {
+      aArr[i]   = new double [Nz];
+      vArr[i]   = new double [Nz];
+      dArr[i]   = new double [Nz];
+      DpdArr[i] = new double [Nz];
+      DzdArr[i] = new double [Nz];
+    }
+  } else {
+    Np = 0; ip = 0.0; fp = 0.0; dp = 1.0;
+    Nz = 0; iz = 0.0; fz = 0.0; dz = 1.0;
+    TFO = 0.0;
+    aArr   = NULL;
+    vArr   = NULL;
+    dArr   = NULL;
+    DpdArr = NULL;
+    DzdArr = NULL;
+  }
+  delete HSFile;
+
+  // Read FOHSa.txt file and fill array
+  FName = "FOHSa.txt";
+  HSFile = new ifstream;
+  HSFile->open(FName);
+  if(HSFile->is_open()) {
+    HSFile->seekg (0, ios::beg);
+    for(i=0;i<Np;i++) {
+      for(j=0;j<Nz;j++) {
+        HSFile->getline(buff,100);
+        aArr[i][j]=atof(buff);
+      }
+    }
+    HSFile->close();
+  }
+  delete HSFile;
+
+  // Read FOHSv.txt file and fill array
+  FName = "FOHSv.txt";
+  HSFile = new ifstream;
+  HSFile->open(FName);
+  if(HSFile->is_open()) {
+    HSFile->seekg (0, ios::beg);
+    for(i=0;i<Np;i++) {
+      for(j=0;j<Nz;j++) {
+        HSFile->getline(buff,100);
+        vArr[i][j]=atof(buff);
+      }
+    }
+    HSFile->close();
+  }
+  delete HSFile;
+
+  // Read FOHSd.txt file and fill array
+  FName = "FOHSd.txt";
+  HSFile = new ifstream;
+  HSFile->open(FName);
+  if(HSFile->is_open()) {
+    HSFile->seekg (0, ios::beg);
+    for(i=0;i<Np;i++) {
+      for(j=0;j<Nz;j++) {
+        HSFile->getline(buff,100);
+        dArr[i][j]=atof(buff);
+      }
+    }
+    HSFile->close();
+  }
+  delete HSFile;
+
+  // Read FOHSDpd.txt file and fill array
+  FName = "FOHSDpd.txt";
+  HSFile = new ifstream;
+  HSFile->open(FName);
+  if(HSFile->is_open()) {
+    HSFile->seekg (0, ios::beg);
+    for(i=0;i<Np;i++) {
+      for(j=0;j<Nz;j++) {
+        HSFile->getline(buff,100);
+        DpdArr[i][j]=atof(buff);
+      }
+    }
+    HSFile->close();
+  }
+  delete HSFile;
+
+  // Read FOHSDzd.txt file and fill array
+  FName = "FOHSDzd.txt";
+  HSFile = new ifstream;
+  HSFile->open(FName);
+  if(HSFile->is_open()) {
+    HSFile->seekg (0, ios::beg);
+    for(i=0;i<Np;i++) {
+      for(j=0;j<Nz;j++) {
+        HSFile->getline(buff,100);
+        DzdArr[i][j]=atof(buff);
+      }
+    }
+    HSFile->close();
+  }
+  delete HSFile;
+}
+
+Hypersurface::~Hypersurface(void) {
+  for(i=0;i<Np;i++) {
+    delete[] aArr[i];
+    delete[] vArr[i];
+    delete[] dArr[i];
+    delete[] DpdArr[i];
+    delete[] DzdArr[i];
+  }
+  delete[] aArr;
+  delete[] vArr;
+  delete[] dArr;
+  delete[] DpdArr;
+  delete[] DzdArr;
+}
+
+double Hypersurface::fahs(double p, double z) {
+  i=int(p/dp);
+  j=int(z/dz);
+  if(i>=Np-1) i=Np-2;
+  if(j>=Nz-1) j=Nz-2;
+  return (aArr[i][j]   * (i+1-p/dp) + aArr[i+1][j]   * (p/dp-i)) * (j+1-z/dz) +
+         (aArr[i][j+1] * (i+1-p/dp) + aArr[i+1][j+1] * (p/dp-i)) * (z/dz-j);
+}
+
+double Hypersurface::fvhs(double p, double z) {
+  i=int(p/dp);
+  j=int(z/dz);
+  if(i>=Np-1) i=Np-2;
+  if(j>=Nz-1) j=Nz-2;
+  return (vArr[i][j]   * (i+1-p/dp) + vArr[i+1][j]   * (p/dp-i)) * (j+1-z/dz) +
+         (vArr[i][j+1] * (i+1-p/dp) + vArr[i+1][j+1] * (p/dp-i)) * (z/dz-j);
+}
+
+double Hypersurface::fdhs(double p, double z) {
+  i=int(p/dp);
+  j=int(z/dz);
+  if(i>=Np-1) i=Np-2;
+  if(j>=Nz-1) j=Nz-2;
+  return (dArr[i][j]   * (i+1-p/dp) + dArr[i+1][j]   * (p/dp-i)) * (j+1-z/dz) +
+         (dArr[i][j+1] * (i+1-p/dp) + dArr[i+1][j+1] * (p/dp-i)) * (z/dz-j);
+}
+
+double Hypersurface::fDpdhs(double p, double z) {
+  i=int(p/dp);
+  j=int(z/dz);
+  if(i>=Np-1) i=Np-2;
+  if(j>=Nz-1) j=Nz-2;
+  return (DpdArr[i][j]   * (i+1-p/dp) + DpdArr[i+1][j]   * (p/dp-i)) * (j+1-z/dz) +
+         (DpdArr[i][j+1] * (i+1-p/dp) + DpdArr[i+1][j+1] * (p/dp-i)) * (z/dz-j);
+}
+
+double Hypersurface::fDzdhs(double p, double z) {
+  i=int(p/dp);
+  j=int(z/dz);
+  if(i>=Np-1) i=Np-2;
+  if(j>=Nz-1) j=Nz-2;
+  return (DzdArr[i][j]   * (i+1-p/dp) + DzdArr[i+1][j]   * (p/dp-i)) * (j+1-z/dz) +
+         (DzdArr[i][j+1] * (i+1-p/dp) + DzdArr[i+1][j+1] * (p/dp-i)) * (z/dz-j);
+}
diff --git a/TTherminator/Therminator/Hypersurface.h b/TTherminator/Therminator/Hypersurface.h
new file mode 100644 (file)
index 0000000..5981f6b
--- /dev/null
@@ -0,0 +1,58 @@
+/***************************************************************************
+ *                             L H Y Q U I D                               *
+ *              reLativistic HYdrodynamics of Quark glUon fluID            *
+ *                                                                         *
+ *                         T H E R M I N A T O R                           *
+ *                      THERMal heavy-IoN generATOR                        *
+ *                                                                         *
+ *                              INTERFACE                                  *
+ *                                                                         *
+ * Author of LHYQUID THERMINATOR INTERFACE code:                           *
+ *   Miko´┐Żaj Chojnacki, e-mail: Mikolaj.Chojnacki@ifj.edu.pl               *
+ *                                                                         *
+ * About the interface:                                                    *
+ *   Code in the files "Hypersurface.h" and "Hypersurface.cxx" enables the *
+ * usage of a paramatrized hypersurface from LHYQUID (written in           *
+ * Mathematica).                                                           *
+ * Set the option FreezeOutModel = Lhyquid in the file "therminator.in".   *
+ * Modifications in the original files (Integrator.cxx, Integrator.h,      *
+ * therm_events, Makefile, therminator.in are marked by "MCH" comment      *
+ *                                                                         *
+ ***************************************************************************/
+
+
+#ifndef _MCH_HYPERSURFACE_
+#define _MCH_HYPERSURFACE_
+
+#include <stdio.h>
+#include <iostream>
+#include <fstream>
+
+using namespace std;
+
+class Hypersurface {
+  public:
+    // constructor
+    Hypersurface(void);
+    // destructor
+    ~Hypersurface(void);
+    // variables
+    double TFO;
+    // functions
+    double fahs  (double p, double z);
+    double fvhs  (double p, double z);
+    double fdhs  (double p, double z);
+    double fDpdhs(double p, double z);
+    double fDzdhs(double p, double z);
+  private:
+    unsigned int i, j;
+    unsigned int Np, Nz;
+    double   ip, fp, dp, iz, fz, dz;
+    double   **aArr, **vArr, **dArr, **DpdArr, **DzdArr;
+    char     *FName;
+    char     buff[100];
+    ifstream *HSFile;
+};
+
+
+#endif /* _MCH_HYPERSURFACE_ */
diff --git a/TTherminator/Therminator/Integrator.cxx b/TTherminator/Therminator/Integrator.cxx
new file mode 100644 (file)
index 0000000..50edc25
--- /dev/null
@@ -0,0 +1,611 @@
+/******************************************************************************
+ *                      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 <fstream>
+#include "Integrator.h"
+#include "ParticleType.h"
+#include "Particle.h"
+#include <TMath.h>
+#include <TFile.h>
+
+extern ReadPar *sRPInstance;
+extern int      sTables;
+extern int      sModel;
+
+Integrator::Integrator(int aNpart)
+{
+  kFmToGev  = 0.197326960277;                          /*MCH updated: kFmToGev  = 0.197;*/
+  ReadParameters();
+
+  PRINT_MESSAGE("Hash for these parameters is: " << ParameterHash());
+  
+  mNPart = aNpart;
+  kTwoPi2 = TMath::Pi()*TMath::Pi()*2*2;               /*MCH*/
+  kTwoPi3 = TMath::Pi()*TMath::Pi()*TMath::Pi()*2*2*2;
+  mRandom = new TRandom2();
+  //  mRandom->SetSeed2(41321, 8457);
+  mRandom->SetSeed(41321);
+
+  mFOHS = new Hypersurface();                          /*MCH*/
+}
+
+double Integrator::CalcBE(double aX)
+{
+  return 1.0/(TMath::Exp(aX)-1);
+}
+
+double Integrator::CalcFD(double aX)
+{
+  return 1.0/(TMath::Exp(aX)+1);
+}
+
+double Integrator::Calka(double aMass, double aMiu, 
+                        double *aRap, double *aPt, double *aPhiP, 
+                        double *aAlfaP, double *aRho, double *aPhiS, double *aTime, double aSpin,
+                        int aKeepPos)
+{
+  // Generate momentum components
+  (*aRap) = mRandom->Rndm() * mRapRange - mRapRange/2.0;
+  double tZet   = mRandom->Rndm();
+  (*aPt) = tZet/(1-tZet);
+  (*aPhiP) = mRandom->Rndm() * TMath::Pi() * 2.0;
+  
+  // Generate position components
+  if (!aKeepPos) {
+    (*aAlfaP) = mRandom->Rndm()*mAlfaRange - mAlfaRange/2.0;
+    (*aRho) = mRandom->Rndm()*mRhoMax;
+    (*aPhiS) = mRandom->Rndm()*TMath::Pi()*2;
+    if (sModel == 8) {
+      (*aTime) = mTau - mBWDelay*TMath::Log(mRandom->Rndm());
+    }
+  }
+  
+  double tFpod = 0.0;
+  if ((sModel == 0) || (sModel == 3)) { // Single FreezeOut
+    double tPU = (TMath::Sqrt(aMass*aMass+(*aPt)*(*aPt))*
+                cosh((*aAlfaP)-(*aRap))*
+                TMath::Sqrt(1+(((*aRho)*(*aRho))/(mTau*mTau)))
+                -(*aPt)*cos((*aPhiS)-(*aPhiP))*(*aRho)/mTau);
+    if (fabs(floor(aSpin) - aSpin)< 0.01) 
+      tFpod = 1/((1-tZet)*(1-tZet))*mTau*(*aPt)*(*aRho)*tPU*CalcBE((tPU- aMiu)/mTemp)/(kTwoPi3);
+    else
+      tFpod = 1/((1-tZet)*(1-tZet))*mTau*(*aPt)*(*aRho)*tPU*CalcFD((tPU- aMiu)/mTemp)/(kTwoPi3);
+  }
+/*MCH begin*/
+  else if (sModel == 10) { // Lhyquid3D
+    double tZeta  = mRandom->Rndm()*0.5*TMath::Pi();           // angle in rho-t plane
+    double taHS   = mFOHS->fahs  ((*aPhiS),tZeta);             // velocity angle; 0.0=radial flow
+    double tvHS   = mFOHS->fvhs  ((*aPhiS),tZeta);             // velocity
+    double tdHS   = mFOHS->fdhs  ((*aPhiS),tZeta)/kFmToGev;    // distance in rho-t plane           [GeV^-1]
+    double tDpdHS = mFOHS->fDpdhs((*aPhiS),tZeta)/kFmToGev;    // distance derivative over (*aPhiS) [GeV^-1]
+    double tDzdHS = mFOHS->fDzdhs((*aPhiS),tZeta)/kFmToGev;    // distance derivative over tZeta    [GeV^-1]
+    double tTemp  = mFOHS->TFO;                                        // freeze-out temparature
+    double tMt    = TMath::Hypot(aMass, (*aPt));               // transverse mass
+    (*aRho)       = tdHS*cos(tZeta);                           // rho
+    double tdPt   = 1.0/((1-tZet)*(1-tZet));                   // dPt
+    (*aTime)      = tdHS*sin(tZeta)*cosh(*aAlfaP);             // t
+
+    double tPU    = 1.0/sqrt(1.0-tvHS*tvHS)*(tMt*cosh((*aRap)-(*aAlfaP))-(*aPt)*tvHS*cos((*aPhiS)-(*aPhiP)+taHS));
+    double tFC    = tdHS*tdHS*sin(tZeta)*(
+                     tdHS  *cos(tZeta)*( tMt*sin(tZeta)*cosh((*aRap)-(*aAlfaP))+(*aPt)*cos(tZeta)*cos((*aPhiS)-(*aPhiP)))+
+                     tDzdHS*cos(tZeta)*(-tMt*cos(tZeta)*cosh((*aRap)-(*aAlfaP))+(*aPt)*sin(tZeta)*cos((*aPhiS)-(*aPhiP)))+
+                     tDpdHS*(*aPt)*sin((*aPhiS)-(*aPhiP))
+                   );
+    if(tFC < 0.0) tFC = 0.0;
+    if (fabs(floor(aSpin)-aSpin) < 0.01)
+      tFpod =1.0/kTwoPi3 * tdPt*(*aPt) * tFC * CalcBE((tPU-aMiu)/tTemp);
+    else
+      tFpod =1.0/kTwoPi3 * tdPt*(*aPt) * tFC * CalcFD((tPU-aMiu)/tTemp);
+  }
+  else if (sModel == 11) { // Lhyquid2D
+    (*aAlfaP)     = (*aRap);                                   // dirac delta (Y-ap)
+    double tZeta  = mRandom->Rndm()*0.5*TMath::Pi();           // angle in rho-t plane
+    double taHS   = mFOHS->fahs  ((*aPhiS),tZeta);             // velocity angle; 0.0=radial flow
+    double tvHS   = mFOHS->fvhs  ((*aPhiS),tZeta);             // velocity
+    double tdHS   = mFOHS->fdhs  ((*aPhiS),tZeta)/kFmToGev;    // distance in rho-t plane           [GeV^-1]
+    double tDpdHS = mFOHS->fDpdhs((*aPhiS),tZeta)/kFmToGev;    // distance derivative over (*aPhiS) [GeV^-1]
+    double tDzdHS = mFOHS->fDzdhs((*aPhiS),tZeta)/kFmToGev;    // distance derivative over tZeta    [GeV^-1]
+    double tTemp  = mFOHS->TFO;                                        // freeze-out temparature
+    double tMt    = TMath::Hypot(aMass, (*aPt));               // transverse mass
+    (*aRho)       = tdHS*cos(tZeta);                           // rho
+    double tdPt   = 1.0/((1-tZet)*(1-tZet));                   // dPt
+    (*aTime)      = tdHS*sin(tZeta)*cosh(*aAlfaP);             // t
+
+    double tPU    = 1.0/sqrt(1.0-tvHS*tvHS)*(tMt-(*aPt)*tvHS*cos((*aPhiS)-(*aPhiP)+taHS));
+    double tFC    = tdHS*(
+                     tdHS  *cos(tZeta)*( (*aPt)/tMt*cos(tZeta)*cos((*aPhiS)-(*aPhiP))+sin(tZeta) )+
+                     tDzdHS*cos(tZeta)*( (*aPt)/tMt*sin(tZeta)*cos((*aPhiS)-(*aPhiP))-cos(tZeta) )+
+                     tDpdHS*             (*aPt)/tMt*           sin((*aPhiS)-(*aPhiP))
+                   );
+    if(tFC < 0.0) tFC = 0.0;
+    if (fabs(floor(aSpin)-aSpin) < 0.01)
+      tFpod =1.0/kTwoPi2 * tdPt*(*aPt) * tFC * CalcBE((tPU-aMiu)/tTemp);
+    else
+      tFpod =1.0/kTwoPi2 * tdPt*(*aPt) * tFC * CalcFD((tPU-aMiu)/tTemp);
+  }
+/*MCH end*/  
+  else if (sModel == 2) { // Blast-Wave with Vt
+    double tMt  = TMath::Hypot(aMass, (*aPt));
+    double tPre = (mTau + mBWA*(*aRho))*(*aPt)*(*aRho)/kTwoPi3;
+    double tCHay = cosh((*aAlfaP)-(*aRap));
+    double tCSsp = cos((*aPhiS) - (*aPhiP));
+    double tBra = tMt*tCHay - mBWA*(*aPt)*tCSsp;
+    double tGamma = 1.0/sqrt(1-mBWVt*mBWVt);
+    double tPU  = tMt*tGamma*tCHay-(*aPt)*mBWVt*tGamma*tCSsp - aMiu;
+    if (fabs(floor(aSpin) - aSpin)< 0.01) 
+      tFpod = 1/((1-tZet)*(1-tZet))*tPre*tBra*CalcBE(tPU/mTemp);
+    else
+      tFpod = 1/((1-tZet)*(1-tZet))*tPre*tBra*CalcFD(tPU/mTemp);
+  }
+  else if (sModel == 6) { // Blast-Wave with Vt
+    double tMt  = TMath::Hypot(aMass, (*aPt));
+    double tPre = (mTau + mBWA*(*aRho))*(*aPt)*(*aRho)/kTwoPi3;
+    double tCHay = cosh((*aAlfaP)-(*aRap));
+    double tCSsp = cos((*aPhiS) - (*aPhiP));
+    double tBra = tMt*tCHay - mBWA*(*aPt)*tCSsp;
+    double tGamma = 1.0/sqrt(1-mBWVt*mBWVt);
+    double tPU  = tMt*tGamma*tCHay-(*aPt)*mBWVt*tGamma*tCSsp - aMiu;
+    if (fabs(floor(aSpin) - aSpin)< 0.01) 
+      tFpod = 1/((1-tZet)*(1-tZet))*tPre*tBra*CalcBE(tPU/mTemp);
+    else
+      tFpod = 1/((1-tZet)*(1-tZet))*tPre*tBra*CalcFD(tPU/mTemp);
+  }
+  else if (sModel == 4) { // Blast-Wave with linear Vt profile
+    double tMt  = TMath::Hypot(aMass, (*aPt));
+    double tPre = (mTau + mBWA*(*aRho))*(*aPt)*(*aRho)/kTwoPi3;
+    double tCHay = cosh((*aAlfaP)-(*aRap));
+    double tCSsp = cos((*aPhiS) - (*aPhiP));
+    double tBra = tMt*tCHay - mBWA*(*aPt)*tCSsp;
+    double tVt = ((*aRho)/mRhoMax)/(mBWVt+((*aRho)/mRhoMax));
+    double tGamma = 1.0/sqrt(1-tVt*tVt);
+    double tPU  = tMt*tGamma*tCHay-(*aPt)*tVt*tGamma*tCSsp - aMiu;
+    if (fabs(floor(aSpin) - aSpin)< 0.01) 
+      tFpod = 1/((1-tZet)*(1-tZet))*tPre*tBra*CalcBE(tPU/mTemp);
+    else
+      tFpod = 1/((1-tZet)*(1-tZet))*tPre*tBra*CalcFD(tPU/mTemp);
+  }
+  else if (sModel == 7) { 
+    // Blast-Wave with linear Vt profile
+    // and delay in particle emission point - formation time
+    double tMt  = TMath::Hypot(aMass, (*aPt));
+    double tPre = (mTau + mBWA*(*aRho))*(*aPt)*(*aRho)/kTwoPi3;
+    double tCHay = cosh((*aAlfaP)-(*aRap));
+    double tCSsp = cos((*aPhiS) - (*aPhiP));
+    double tBra = tMt*tCHay - mBWA*(*aPt)*tCSsp;
+    double tVt = ((*aRho)/mRhoMax)/(mBWVt+((*aRho)/mRhoMax));
+    double tGamma = 1.0/sqrt(1-tVt*tVt);
+    double tPU  = tMt*tGamma*tCHay-(*aPt)*tVt*tGamma*tCSsp - aMiu;
+    if (fabs(floor(aSpin) - aSpin)< 0.01) 
+      tFpod = 1/((1-tZet)*(1-tZet))*tPre*tBra*CalcBE(tPU/mTemp);
+    else
+      tFpod = 1/((1-tZet)*(1-tZet))*tPre*tBra*CalcFD(tPU/mTemp);
+  }
+  else if (sModel == 8) { 
+    // Blast-Wave with linear Vt profile
+    // and delay in emission time - exponential decay
+    double tTau   = (*aTime);
+    double tMt    = TMath::Hypot(aMass, (*aPt));
+    double tPre   = (tTau + mBWA*(*aRho))*(*aPt)*(*aRho)/kTwoPi3;
+    double tCHay  = cosh((*aAlfaP)-(*aRap));
+    double tCSsp  = cos((*aPhiS) - (*aPhiP));
+    double tBra   = tMt*tCHay - mBWA*(*aPt)*tCSsp;
+    double tVt    = ((*aRho)/mRhoMax)/(mBWVt+((*aRho)/mRhoMax));
+    double tGamma = 1.0/sqrt(1-tVt*tVt);
+    double tPU    = tMt*tGamma*tCHay-(*aPt)*tVt*tGamma*tCSsp - aMiu;
+    if (fabs(floor(aSpin) - aSpin)< 0.01) 
+      tFpod = 1/((1-tZet)*(1-tZet))*tPre*tBra*CalcBE(tPU/mTemp);
+    else
+      tFpod = 1/((1-tZet)*(1-tZet))*tPre*tBra*CalcFD(tPU/mTemp);
+  }
+/* -=[ Testing with Mathematica ]=-*/
+/*if ((sModel == 10)||(sModel == 11)) {
+  cout.precision(20);
+  cout << endl;
+  cout << "{phi  = " << (*aPhiS)  << ","  << endl;
+  cout << " zeta = " <<   tZeta   << ","  << endl;
+  cout << " z    = " <<   tZet    << ","  << endl;
+  cout << " phip = " << (*aPhiP)  << ","  << endl;
+  cout << " ap   = " << (*aAlfaP) << ","  << endl;
+  cout << " Y    = " << (*aRap)   << ","  << endl;
+  cout << " m    = " <<   aMass   << ","  << endl;
+  cout << " s    = " <<   aSpin   << ","  << endl;
+  cout << " mu   = " <<   aMiu    << "};" << endl;
+  cout << endl;
+  cout.precision(6);
+  cout << "TF0  = " <<   tTemp  << endl;
+  cout << "t    = " << (*aTime) << endl;
+  cout << "rho  = " << (*aRho)  << endl;
+  cout << "pT   = " << (*aPt)   << endl;
+  cout << "dpT  = " <<   tdPt   << endl;
+  cout << "mT   = " <<   tMt    << endl;
+  cout << endl;
+  cout << "aHS     = " << taHS     << endl;
+  cout << "vHS     = " << tvHS     << endl;
+  cout << "dHS     = " << tdHS     << endl;
+  cout << "DpdHS   = " << tDpdHS   << endl;
+  cout << "DzdHS   = " << tDzdHS   << endl;
+  cout << "PU      = " << tPU      << endl;
+  cout << "FC      = " << tFC      << endl;
+  cout << "part1   = " << (tPU - aMiu)/tTemp << endl;
+  cout << "part2BE = " << CalcBE( (tPU - aMiu)/tTemp ) << endl;
+  cout << "part2FD = " << CalcFD( (tPU - aMiu)/tTemp ) << endl;
+  cout << "Fpod    = " << tFpod    << endl;
+  getchar();
+}*/
+
+  return tFpod;
+}
+
+double Integrator::GetMiu(double aIzo, double aBar, double aStr)
+{
+  return (mMiu_i*aIzo+mMiu_s*aStr+mMiu_b*aBar);
+}
+
+double Integrator::CalcFunPodCalk(double aMass, double aIzo, double aBar, double aStr, double aSpin)
+{
+  //  int proc = mNPart/100;
+  double tFpod;
+  double tMax = 0.0;
+  double tMiu = GetMiu(aIzo, aBar, aStr);
+  double tRap, tPhiP, tAlfaP, tRho, tPhiS, tTime; 
+  double tPt;
+
+  for (int tIpart=0; tIpart<mNPart; tIpart++)
+    {
+      tFpod = Calka(aMass,tMiu,&tRap,&tPt,&tPhiP,&tAlfaP,&tRho,&tPhiS,&tTime,aSpin);
+      if (tFpod>tMax) tMax = tFpod;
+    }
+  return tMax;
+}
+
+double Integrator::Integrate(double aMass, double aIzo, double aBar, double aStr, double aSpin)
+{
+  double tFpod=0.0;
+  double tFtest=0.0;
+  double tMiu = GetMiu(aIzo, aBar, aStr);
+  double tRap, tPhiP, tAlfaP, tRho, tPhiS, tTime; 
+  double tPt;
+
+  for (int tIpart=0; tIpart<mNPart; tIpart++)
+    {
+      tFpod = Calka(aMass,tMiu,&tRap,&tPt,&tPhiP,&tAlfaP,&tRho,&tPhiS,&tTime,aSpin);
+
+      tFtest += tFpod;
+    }
+
+  double tCalka;
+  if (sModel == 5)
+    tCalka = mRapRange*(1.0)*mAlfaRange*4*TMath::Pi()*TMath::Pi()*(mTauf-mTau0)* tFtest / (mNPart*1.0);  
+/*MCH begin*/
+  else if (sModel == 10)
+    tCalka = mRapRange*mAlfaRange*(1.0)*2.0*TMath::Pi()*2.0*TMath::Pi()*0.5*TMath::Pi()*tFtest / (mNPart*1.0);
+  else if (sModel == 11)
+    tCalka = mRapRange*(1.0)*2.0*TMath::Pi()*2.0*TMath::Pi()*0.5*TMath::Pi()*tFtest / (mNPart*1.0);
+/*MCH end*/    
+  else
+    tCalka = mRapRange*(1.0)*mAlfaRange*4*TMath::Pi()*TMath::Pi()*mRhoMax* tFtest / (mNPart*1.0);  
+  return tCalka;
+}
+
+void 
+Integrator::Generate(ParticleType *aPartType, int aPartCount, Particle ***oParticles)
+{
+  int tIter = 0;
+  double tFpod;
+  double tFtest;
+  double tMiu = GetMiu(1.0*aPartType->GetI3(), 1.0*aPartType->GetBarionN(), 1.0*aPartType->GetStrangeness());
+  double tFMax;
+  double tSpin = aPartType->GetSpin();
+  double tRap, tPhiP, tAlfaP, tRho, tPhiS, tTime; 
+  double tPt;
+
+  PRINT_DEBUG_3("Gen for " << (aPartType->GetName()) << " i3 b s " << 1.0*aPartType->GetI3() << " " << 1.0*aPartType->GetBarionN() << " " << 1.0*aPartType->GetStrangeness() << " " << tMiu);
+  
+  tFMax = aPartType->GetFMax();
+
+  (*oParticles) = (Particle **) malloc(sizeof(Particle *) * aPartCount);
+  Particle *tBuf;
+  
+  while (tIter<aPartCount)
+    {
+      tFpod = Calka(aPartType->GetMass(),tMiu,&tRap,&tPt,&tPhiP,&tAlfaP,&tRho,&tPhiS,&tTime,tSpin);
+      tFtest = mRandom->Rndm()*tFMax;
+      if (tFtest<tFpod)
+       {
+         if ((sModel == 0) || (sModel == 3)) // Single freeze-out
+           tBuf = new Particle(tRap, tPt, tPhiP, tAlfaP, tRho, tPhiS, TMath::Hypot(mTau, tRho), aPartType);
+/*MCH begin*/
+         else if ((sModel == 10) || (sModel == 11)) {
+           tBuf = new Particle(tRap, tPt, tPhiP, tAlfaP, tRho, tPhiS, tTime/cosh(tAlfaP), aPartType);
+         }
+/*MCH end*/
+         else if ((sModel == 1) || (sModel == 2) || (sModel == 4)) { // Blast-wave
+           double tTau = mTau + mBWA * tRho;
+           tBuf = new Particle(tRap, tPt, tPhiP, tAlfaP, tRho, tPhiS, tTau, aPartType);
+         }
+         else if (sModel == 7) { // Blast-wave
+           double tTau = mTau + mBWA * tRho;
+           double px = tPt*TMath::Cos(tPhiP);
+           double py = tPt*TMath::Sin(tPhiP);
+           double tMt = TMath::Hypot(aPartType->GetMass(),tPt);
+           double pz = tMt*TMath::SinH(tRap);
+  
+           double rx = tRho*TMath::Cos(tPhiS);
+           double ry = tRho*TMath::Sin(tPhiS);
+
+           double rz = tTau*TMath::SinH(tAlfaP);
+           double rt = tTau*TMath::CosH(tAlfaP);
+           double dt = -mBWDelay * TMath::Log(mRandom->Rndm());
+           rt += dt;
+           double en = sqrt(tMt*tMt + pz*pz);
+           rx += dt * px/en;
+           ry += dt * py/en;
+           rz += dt * pz/en;
+
+           tBuf = new Particle(aPartType, px, py, pz, rx, ry, rz, rt);
+         }
+         else if (sModel == 8) { // Blast-wave
+           double tTau = tTime + mBWA * tRho;
+           tBuf = new Particle(tRap, tPt, tPhiP, tAlfaP, tRho, tPhiS, tTau, aPartType);
+         }
+         else if (sModel == 6) { // Blast-wave
+           double tTau = mTau + mBWA * tRho;
+           double px = tPt*TMath::Cos(tPhiP);
+           double py = tPt*TMath::Sin(tPhiP);
+           double tMt = TMath::Hypot(aPartType->GetMass(),tPt);
+           double pz = tMt*TMath::SinH(tRap);
+  
+           double rx = tRho*TMath::Cos(tPhiS);
+           double ry = tRho*TMath::Sin(tPhiS);
+
+           double rz = tTau*TMath::SinH(tAlfaP);
+           double rt = tTau*TMath::CosH(tAlfaP);
+           rt += -mBWDelay * TMath::Log(mRandom->Rndm());
+
+           tBuf = new Particle(aPartType, px, py, pz, rx, ry, rz, rt);
+         }
+         else if (sModel == 5) {
+           double tTau = sqrt(tTime*tTime - tAlfaP*tAlfaP);
+           double tAlfa = 0.5*log((tTime+tAlfaP)/(tTime-tAlfaP));
+           tBuf = new Particle(tRap, tPt, tPhiP, tAlfa, tRho, tPhiS, tTau, aPartType);
+         }
+         (*oParticles)[tIter] = tBuf;
+         tIter++;
+       }
+    }
+}
+
+void  
+Integrator::Randomize()
+{
+  TDatime tDat;
+
+  //  mRandom->SetSeed2(tDat.Get(), (tDat.Get() % 11) * 7 + (tDat.Get() / 7));
+  mRandom->SetSeed(tDat.Get());
+}
+
+void  
+Integrator::ReadParameters()
+{
+  STR tModel;
+  STR tTable;
+
+  // First read the model parameters
+  try {
+    
+    mRhoMax = atof(sRPInstance->getPar("RhoMax").Data()) / kFmToGev;
+    mTau    = atof(sRPInstance->getPar("Tau").Data()) / kFmToGev;
+    mTemp   = atof(sRPInstance->getPar("Temperature").Data());
+    mMiu_i  = atof(sRPInstance->getPar("MiuI").Data());
+    mMiu_s  = atof(sRPInstance->getPar("MiuS").Data());
+    mMiu_b  = atof(sRPInstance->getPar("MiuB").Data());
+  }
+  catch (STR tError) {
+    PRINT_DEBUG_1("Integrator::ReadParameters - Caught exception " << tError);
+    PRINT_MESSAGE("Did not find one of the neccessary model parameters.");
+    exit(0);
+  }
+
+  // Read additional parameters for BW
+  if ((sModel == 1) || (sModel == 2) || (sModel == 4) || (sModel == 7) || (sModel == 8)) {
+    try {
+      mBWA  = atof(sRPInstance->getPar("BWA").Data());
+    }
+    catch (STR tError) {
+      // Using default value of 0
+      mBWA = 0.0;
+    }
+  }
+  
+  // Read additional parameters for MBWVt
+  if ((sModel == 2) || (sModel == 4) || (sModel == 7) || (sModel == 8)) {
+    try {
+      mBWVt  = atof(sRPInstance->getPar("BWVt").Data());
+    }
+    catch (STR tError) {
+      PRINT_DEBUG_1("Integrator::ReadParameters (BW Vt part) - Caught exception " << tError);
+      PRINT_MESSAGE("Did not find one of the neccessary model parameters.");
+      exit(0);
+    }
+  }
+
+  // Read additional parameters for MBWVt
+  if ((sModel ==6) || (sModel == 7) || (sModel == 8))
+    try {
+      mBWDelay  = atof(sRPInstance->getPar("BWDelay").Data()) / kFmToGev;
+    }
+    catch (STR tError) {
+      PRINT_DEBUG_1("Integrator::ReadParameters (BW Vt Delay part) - Caught exception " << tError);
+      PRINT_MESSAGE("Did not find one of the neccessary model parameters.");
+      exit(0);
+    }
+  
+  
+  // Then the integration range
+  try {
+    mAlfaRange = atof(sRPInstance->getPar("AlphaRange").Data());
+    mRapRange  = atof(sRPInstance->getPar("RapidityRange").Data());
+  }
+  catch (STR tError) {
+    PRINT_DEBUG_1("Integrator::ReadParameters - Caught exception " << tError);
+    PRINT_MESSAGE("Did not find one of the neccessary integration ranges.");
+    exit(0);
+  }
+
+  // Read hydro-specific parameters
+  if (sModel == 5) {
+    try {
+      mTauf   = atof(sRPInstance->getPar("TauF").Data());
+      mTau0   = atof(sRPInstance->getPar("Tau0").Data());
+      mLambda = atof(sRPInstance->getPar("Lambda").Data());
+      mBN     = atof(sRPInstance->getPar("BN").Data());
+      mAlfa   = atof(sRPInstance->getPar("Alfa").Data());
+    }
+    catch (STR tError) {
+      PRINT_DEBUG_1("Integrator::ReadParameters - Caught exception " << tError);
+      PRINT_MESSAGE("Did not find one of the neccessary hydro parameters.");
+      exit(0);
+    }
+  }
+}
+
+char *  
+Integrator::ParameterHash()
+{
+  char *tBuf;
+  
+  tBuf = (char *) malloc(sizeof(char *) * (15 * 3 + 3));
+  
+  if (sModel == 0) {
+    sprintf(tBuf, "%1i%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f", sTables, mRhoMax*100, mTau*100, mTemp*1000, mMiu_i*10000, mMiu_s*100000, mMiu_b*10000, mAlfaRange*100, mRapRange*100);
+  }
+/*MCH begin*/
+  else if ((sModel == 10) || (sModel == 11)) {
+    sprintf(tBuf, "%1i%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f", sTables, mRhoMax*100, mTau*100, mTemp*1000, mMiu_i*10000, mMiu_s*100000, mMiu_b*10000, mAlfaRange*100, mRapRange*100);
+  }
+/*MCH end*/
+  else if ((sModel == 2) || (sModel == 4))  {
+    sprintf(tBuf, "%1i%1i%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f", sModel, sTables, mRhoMax*100, mTau*100, mTemp*1000, mMiu_i*10000, mMiu_s*100000, mMiu_b*10000, mBWA*100, mBWVt*1000, mAlfaRange*100, mRapRange*100);
+  }
+  else if ((sModel == 6) || (sModel == 7) || (sModel == 8)) {
+    sprintf(tBuf, "%1i%1i%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f", sModel, sTables, mRhoMax*100, mTau*100, mTemp*1000, mMiu_i*10000, mMiu_s*100000, mMiu_b*10000, mBWA*100, mBWVt*1000, mAlfaRange*100, mRapRange*100, mBWDelay);
+  }
+  else if (sModel == 5)  {
+    sprintf(tBuf, "%1i%1i%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f", sModel, sTables, mRhoMax*100, mTemp*1000, mMiu_i*10000, mMiu_s*100000, mMiu_b*10000, mTau0*10, mTauf*10, mAlfa*1000, mBN*100, mLambda*100, mAlfaRange*100, mRapRange*100);
+  }
+  
+  return tBuf;
+}
+
+void   
+Integrator::ReadMultInteg(ParticleDB *aDB)
+{
+  // Make or read table with propabilities
+  int tK;
+  char *tHash;
+  char  tIntName[100];
+  char  tMultName[100];
+  ifstream *tFileIn = NULL;
+  ofstream *tFileOut = NULL;
+  
+  tHash = ParameterHash();
+  
+  strcpy(tIntName, "fintegrandmax_");
+  strcat(tIntName, tHash);
+  strcat(tIntName, ".txt");
+  
+  tFileIn = new ifstream(tIntName);
+  if ((tFileIn) && (tFileIn->is_open())) {
+    PRINT_MESSAGE("Reading Max Integrand values from " << tIntName);
+
+    char tPart[255];
+    float tRFunMax;
+    std::string tStrPart;
+    
+    while (!tFileIn->eof()) {
+      (*tFileIn) >> tPart >> tRFunMax;
+      
+      tStrPart = tPart;
+      aDB->GetParticleType(tStrPart)->SetFMax(tRFunMax);
+    }
+  }
+  else {
+    float tRFunMax;
+
+    PRINT_MESSAGE("Max Integrand file " << tIntName << " not found. Generating...");
+    
+    tFileOut = new ofstream(tIntName);
+
+    for(tK=0; tK<aDB->GetParticleTypeCount(); tK++)
+      {
+       tRFunMax = CalcFunPodCalk(aDB->GetParticleType(tK)->GetMass(),
+                              aDB->GetParticleType(tK)->GetI3(),
+                              aDB->GetParticleType(tK)->GetBarionN()*1.0,
+                              aDB->GetParticleType(tK)->GetStrangeness()*1.0,
+                              aDB->GetParticleType(tK)->GetSpin());
+       
+       (*tFileOut) << aDB->GetParticleType(tK)->GetName() << " "
+               << tRFunMax <<endl;
+       aDB->GetParticleType(tK)->SetFMax(tRFunMax);
+       PRINT_DEBUG_1("particle "<<aDB->GetParticleType(tK)->GetNumber()<<":"
+                     <<aDB->GetParticleType(tK)->GetName()
+                     <<" done");
+      }
+    tFileOut->close();
+  }
+  
+  // Calculate or read multiplicities
+  strcpy(tMultName, "fmultiplicity_");
+  strcat(tMultName, tHash);
+  strcat(tMultName, ".txt");
+
+  tFileIn = new ifstream(tMultName);
+  if ((tFileIn) && (tFileIn->is_open())) {
+    PRINT_MESSAGE("Reading Multiplicities from " << tMultName);
+  }
+  else {
+    PRINT_MESSAGE("Multiplicities file " << tMultName << " not found. Generating...");
+    
+    tFileOut = new ofstream(tMultName);
+      
+    for (int tPart=0; tPart<aDB->GetParticleTypeCount(); tPart++) {
+      double tMult = Integrate(aDB->GetParticleType(tPart)->GetMass(),
+                              aDB->GetParticleType(tPart)->GetI3(),
+                              aDB->GetParticleType(tPart)->GetBarionN()*1.0,
+                              aDB->GetParticleType(tPart)->GetStrangeness()*1.0,
+                              aDB->GetParticleType(tPart)->GetSpin());
+      (*tFileOut) << (aDB->GetParticleType(tPart)->GetName()) << " " << tMult*TMath::Abs(aDB->GetParticleType(tPart)->GetSpin()*2+1) << endl;
+
+      PRINT_DEBUG_1("particle "<<aDB->GetParticleType(tPart)->GetNumber()<<":"
+                   <<aDB->GetParticleType(tPart)->GetName()
+                   <<" done");
+    }
+    tFileOut->close();
+  }
+      
+  free (tHash);
+}
+
diff --git a/TTherminator/Therminator/Integrator.h b/TTherminator/Therminator/Integrator.h
new file mode 100644 (file)
index 0000000..67dae7c
--- /dev/null
@@ -0,0 +1,102 @@
+/******************************************************************************
+ *                      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.                                                                 *
+ *                                                                            *
+ *****************************************************************************/
+#ifndef _BFPW_INTEGRATOR_
+#define _BFPW_INTEGRATOR_
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <iostream>
+#include <math.h>
+#include <string>
+#include <fstream>
+#include <TRandom2.h>
+#include "THGlobal.h"
+#include "ReadPar.h"
+#include "ParticleDB.h"
+#include "Hypersurface.h"      /*MCH*/
+
+class ParticleType;
+class Particle;
+
+class Integrator
+{
+ public:
+  Integrator(int);
+  ~Integrator(){}
+         
+  double CalcBE(double);
+  double CalcFD(double);
+  double GetMiu(double,double,double);
+
+  double Calka(double,double,double*,double*,double*,double*,double*,double*,double*,double,int keeppos=0);
+  double CalcFunPodCalk(double,double,double,double,double);
+  double Integrate(double,double,double,double,double);
+  void   Generate(ParticleType *aPartType, int aPartCount, Particle*** oParticles);
+  
+  void   Randomize();
+  void   ReadMultInteg(ParticleDB *aDB);
+  char  *ParameterHash();
+
+ private:
+  void   ReadParameters();
+  
+  int    mNPart;
+  double kFmToGev;
+  double mRhoMax;
+  double mTau;
+  double mTemp;
+  double mMiu_i;
+  double mMiu_s;
+  double mMiu_b;
+
+  double mBWA;
+  double mBWP;
+
+  double mBWVt;
+  
+  double mBWDelay;
+
+  double mAlfaRange;
+  double mRapRange;
+  
+  double kTwoPi2;              /*MCH*/
+  double kTwoPi3;
+
+  TRandom2 *mRandom;
+
+  // Hydro parameters
+  double mTauf;
+  double mTau0;
+  double mLambda;
+  double mBN;
+  double mAlfa;
+
+  Hypersurface *mFOHS;         /*MCH*/
+};
+
+#endif
+
diff --git a/TTherminator/Therminator/Makefile b/TTherminator/Therminator/Makefile
new file mode 100644 (file)
index 0000000..460fe16
--- /dev/null
@@ -0,0 +1,42 @@
+CC      = g++
+FF      = g77
+LD      = g++
+
+BINARIES = therm_events therm_tree
+HSOURCES = Parser.cxx ParticleType.cxx DecayChannel.cxx Integrator.cxx ParticleDB.cxx DecayTable.cxx Particle.cxx Event.cxx ParticleDecayer.cxx ReadPar.cxx Hypersurface.cxx
+HEADERS  = $(HSOURCES:.cxx=.h)
+SOURCES  = $(HSOURCES) therm_events.cxx
+OBJECTS  = $(SOURCES:.cxx=.o) 
+
+CPPOPT   = -Wno-deprecated `root-config --cflags`
+COPTN    = -c -O0 -g ${CPPOPT} -D_DEBUG_LEVEL_=1
+COPT     = -c -O0 -g ${CPPOPT} -D_DEBUG_LEVEL_=1
+LLIB     = `root-config --libs` -L$(ROOTSYS)/lib -lm -lgcc -g
+
+all: $(BINARIES)
+
+therm_events: $(OBJECTS)
+       $(LD) $(LLIB) -o $@ $(OBJECTS) 
+
+therm_tree: therm_tree.o
+       $(LD) $(LLIB) -o $@ therm_tree.o
+
+therm_tree.o: therm_tree.C
+       $(CC) $(COPT) -o $@ $^
+
+%.o: %.cxx
+       $(CC) $(COPT) -o $@ $< 
+
+%.d: %.cxx
+       $(SHELL) -ec '$(CC) -MM $(CPPOPT) $< | sed '\''s/\($*\)\.o[ :]*/\1.o $@ : /g'\'' > $@; [ -s $@ ] || rm -f $@'
+
+include $(SOURCES:.cxx=.d)
+
+clean:
+       rm -f *.o therm_events therm_tree therm_events.exe therm_tree.exe *.d 
+
+cleanod:
+       rm -f *.o *.d 
+
+package: $(SOURCES) $(HEADERS) Makefile  therminator.in therm_tree.C ./share/particles.data ./share/decays.data THGlobal.h figure*.C runevents.sh
+       tar zcvf therminator-2_0.tar.gz $^
diff --git a/TTherminator/Therminator/Parser.cxx b/TTherminator/Therminator/Parser.cxx
new file mode 100644 (file)
index 0000000..78d154d
--- /dev/null
@@ -0,0 +1,775 @@
+/******************************************************************************
+ *                      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 <TMath.h>
+#include "THGlobal.h"
+#include "ReadPar.h"
+#include "Parser.h"
+#include "DecayChannel.h"
+#include <sstream>
+
+extern ReadPar *sRPInstance;
+
+const double factorials[7] = {1.0, 1.0, 2.0, 6.0, 24.0, 120.0, 720.0 };
+
+Parser::Parser(ParticleDB *aDB)
+{
+  mDB = aDB;
+  ReadParameters();
+}
+
+Parser::~Parser()
+{
+}
+
+int Parser::check(char *a,char *b)
+{
+  int tIter3=0;
+
+  while(a[tIter3]!='\0') 
+    {
+      if(a[tIter3]!=b[tIter3]) return 0;
+      tIter3++;
+    }
+  return 1;
+}
+
+int Parser::GetParticleCount()
+{
+  return mParticleCount;
+}
+
+char * Parser::GetParticleName(int i)
+{
+  return mNameBuffer[i];
+}
+
+double Parser::GetParticleMass(int i)
+{
+  return mMassBuffer[i];
+}
+
+double Parser::GetParticleGamma(int i)
+{
+  return mGammaBuffer[i];
+}
+
+double Parser::GetParticleSpin(int i)
+{
+  return mSpinBuffer[i];
+}
+
+int Parser::GetParticleBarionN(int i)
+{
+  return mBarionBuffer[i];
+}
+
+int Parser::GetParticleStrangeN(int i)
+{
+  return mStrangeBuffer[i];
+}
+
+double Parser::GetParticleI3(int i)
+{
+  return mI3Buffer[i];
+}
+
+int Parser::GetParticleDecayChannelCount(int i,int j)
+{
+  return mDecayBuffer[i][j];
+}
+
+int Parser::GetParticleNumber(int i)
+{
+  return mTypeCountBuffer[i];
+}
+
+void Parser::ReadInput()
+{
+  int j,tPartIter,l,tIter2, tIter; //variables
+  char str[50];
+  char str1[20];
+  double spin1,spin2,value;
+
+  //////  START OF "TABLES.M" /////////
+  ifstream in("tables.m");
+  if(in)
+    {
+      //START OF HEAD-LINE
+      in.ignore(100,'\n');
+      //END OF HEAD-LINE
+
+      for(tIter2=0;tIter2<50;tIter2++) str[tIter2]='\0';
+
+      tPartIter=0;
+
+      //START OF DATA
+      while(in>>str)
+       {
+         if(*str == 'x') break;
+
+         l=0;
+         //PARTICLE NAME AND MASS
+         if(str[0] == 'm' && str[1] == '[')
+           {
+             for(;;)
+               {
+                 if(str[l+2] == ']') break;
+                 mNameBuffer[tPartIter][l]=str[l+2];   //name
+                 l++;
+               }
+             //                mNameBuffer[tPartIter][l]='\0';
+             in>>str;  // sign "="
+             in>>value;        // mass
+             mMassBuffer[tPartIter]=value;
+             mTypeCountBuffer[tPartIter]=tPartIter;
+             for(tIter2=0;tIter2<50;tIter2++) str[tIter2]='\0';
+           }
+           
+         if(str[0] == 'G' && str[3] == '[')
+           {
+             in>>str;  // sign "="
+             in>>value;        // gamma
+             mGammaBuffer[tPartIter]=value;
+             for(tIter2=0;tIter2<50;tIter2++) str[tIter2]='\0';
+           }
+           
+         // spin
+         if(str[0] == 'J' && str[1] == '[')
+           {
+             in>>str;  // sign "="
+             in>>str;
+             if(str[1] == '/')
+               {
+                 *str1=str[0];spin1=atof(str1);
+                 *str1=str[2];spin2=atof(str1);
+                 mSpinBuffer[tPartIter]=spin1/spin2;
+               }
+             if(str[0] == '-')
+               {
+                 *str1=str[1];spin1=atof(str1);
+                 mSpinBuffer[tPartIter]=-spin1;
+               }
+             if(str[0]!='-' && str[1]!='/') 
+               {
+                 *str1=str[0];
+                 spin1=atof(str1);
+                 mSpinBuffer[tPartIter]=spin1;
+               }
+             tPartIter++;      //next particle
+             for(tIter2=0;tIter2<50;tIter2++) str[tIter2]='\0';
+           }
+       }
+      //END OF DATA
+    }
+  mParticleCount=tPartIter;
+  //////  END OF "TABLES.M" /////////
+
+
+  //////       START OF "i200STAR.m"//////
+
+  double izospin1,izospin2;    //help to calculate izospin, if niecalkowity
+  //    input=fopen("i200STAR.m","r+");
+  ifstream in1("i200STAR.m");
+
+  for(tIter2=0;tIter2<369;tIter2++)
+    for(int jj=0;jj<2;jj++) mDecayBuffer[tIter2][jj]=0; 
+  
+  for(;;)
+    {
+      j=3; //first letter of particle in line in file
+      for(tIter=0;tIter<50;tIter++) str[tIter]='\0';
+      for(tIter=0;tIter<20;tIter++) str1[tIter]='\0';
+      tIter=0;
+       
+      in1.getline(str,200);
+      
+      if(str[0] == 'x') break;
+
+      if(str[0] == 'f' && str[1] == 'i')
+       {       
+         //name
+         for(;;)       
+           {
+             if(str[j] == ',') 
+               {
+                 j++;  //to skip ","
+                 break;
+               }
+             str1[tIter++]=str[j++];
+           }
+         for(tIter=0;tIter<369;tIter++)
+           {       
+             if(check(str1,mNameBuffer[tIter]))
+               {
+                 //barion number
+                 for(tIter2=0;tIter2<20;tIter2++) str1[tIter2]='\0';
+                 if(str[j] == '-')
+                   { 
+                     *str1=str[j+1];
+                     mBarionBuffer[tIter]=atoi(str1);
+                     mBarionBuffer[tIter]=-mBarionBuffer[tIter];
+                   }
+                 if(str[j] != '-') 
+                   {
+                     *str1=str[j];
+                     mBarionBuffer[tIter]=atoi(str1);
+                   }
+                 if(str[j] == '-') j+=3;
+                 else j+=2;
+
+                 //strange number
+                 for(tIter2=0;tIter2<20;tIter2++) str1[tIter2]='\0';
+                 if(str[j] == '-')
+                   {
+                     *str1=str[j+1];
+                     mStrangeBuffer[tIter]=atoi(str1);
+                     mStrangeBuffer[tIter]=-mStrangeBuffer[tIter];
+                   }
+                 if(str[j] != '-') 
+                   {
+                     *str1=str[j];
+                     mStrangeBuffer[tIter]=atoi(str1);
+                   }
+                 if(str[j] == '-') j+=3;
+                 else j+=2;
+
+                 //izospin3 number
+                 for(tIter2=0;tIter2<20;tIter2++) str1[tIter2]='\0';
+                 if(str[j+1] == '/' && str[j+3] == ']') 
+                   {
+                     *str1=str[j];
+                     izospin1=atoi(str1);
+                     *str1=str[j+2];
+                     izospin2=atoi(str1);
+                     mI3Buffer[tIter]=izospin1/izospin2;
+                   }
+                 if(str[j] == '-' && str[j+2] == '/' && str[j+4] == ']') 
+                   {
+                     *str1=str[j+1];
+                     izospin1=atoi(str1);
+                     *str1=str[j+3];
+                     izospin2=atoi(str1);
+                     mI3Buffer[tIter]=-izospin1/izospin2;
+                   }
+                 if(str[j] == '-' && str[j+2] == ']')
+                   {
+                     *str1=str[j+1];
+                     izospin1=atoi(str1);
+                     mI3Buffer[tIter]=-izospin1;
+                   }
+                 if(str[j+1] == ']') 
+                   {
+                     *str1=str[j];
+                     mI3Buffer[tIter]=atof(str1);
+                   }
+                 break;
+               }
+           }
+       }
+       
+      //DECAY CHANNELS
+       
+      tIter=0;
+      //TWO-BODY DECAYS
+      if(str[0] == 's' && str[1] == 'e' && str[2] == '[')
+       {
+         // Reading in the decay channels
+         char *tLBrackert, *tFirstComma, *tSecondComma, *tThirdComma, *tRBracket;
+         char *tFather, *tDaughter1, *tDaughter2, *tBRatio;
+         
+         tLBrackert = strchr(str,'[');
+         tFirstComma = strchr(str,',');
+         tSecondComma = strchr(tFirstComma+1,',');
+         tThirdComma = strchr(tSecondComma+1,',');
+         tRBracket = strchr(tThirdComma,']');
+
+         if (!(tLBrackert && tFirstComma && tSecondComma && tThirdComma && tRBracket))
+           PRINT_DEBUG_1("Malformed line!: " << str);
+         
+         tFather    = strndup(tLBrackert+1,   tFirstComma-tLBrackert-1);
+         tDaughter1 = strndup(tFirstComma+1,  tSecondComma-tFirstComma-1);
+         tDaughter2 = strndup(tSecondComma+1, tThirdComma-tSecondComma-1);
+         tBRatio    = strndup(tThirdComma+1,  tRBracket-tThirdComma-1);
+         
+         // Getting the ratio
+         char *tMiddle, *tRatComponent;
+         double tRatio = 1.0;
+         
+         tMiddle = strchr(tBRatio,'_');
+         if (!tMiddle) tMiddle = strchr(tBRatio,' ');
+         while(tMiddle)
+           {
+             tRatComponent = strndup(tBRatio, tMiddle-tBRatio);
+             if (strchr(tRatComponent,'/'))
+               tRatio *= atof(tRatComponent)/atof(tRatComponent+2);
+             else
+               tRatio *= atof(tRatComponent);
+             tBRatio = tMiddle+1;              
+             tMiddle = strchr(tBRatio,'_');
+             if (!tMiddle) tMiddle = strchr(tBRatio,' ');
+           }
+         if (strchr(tBRatio,'/'))
+           tRatio *= atof(tBRatio)/atof(tBRatio+2);
+         else
+           tRatio *= atof(tBRatio);
+       }
+
+      //THREE-BODY DECAYS
+      j++;     //because se3[]
+      if(str[0] == 's' && str[1] == 'e' && str[2] == '3' && str[3] == '[')
+       {
+         for(;;)       
+           {
+             if(str[j] == ',') 
+               {
+                 j++;  //to skip ","
+                 break;
+               }
+             str1[tIter++]=str[j++];   //name
+           }
+
+         for(tIter=0;tIter<369;tIter++)
+           {       
+             if(check(str1,mNameBuffer[tIter]))
+               {
+                 mDecayBuffer[tIter][1]++;
+                 break;
+               }
+           }
+       }
+
+
+    }
+  //////       END OF "i200STAR.m"//////
+    
+  ParticleType *tPartBuf;
+  int tNum, tPart;
+  
+  for(tPart=0;tPart<mParticleCount;tPart++)
+    {
+      tPartBuf = new ParticleType();
+      tPartBuf->SetName(mNameBuffer[tPart]);
+      tPartBuf->SetMass(mMassBuffer[tPart]);
+      tPartBuf->SetGamma(mGammaBuffer[tPart]);
+      tPartBuf->SetSpin(mSpinBuffer[tPart]);
+      tPartBuf->SetBarionN(mBarionBuffer[tPart]);
+      tPartBuf->SetStrangeness(mStrangeBuffer[tPart]);
+      tPartBuf->SetI3(mI3Buffer[tPart]);
+      tPartBuf->SetDecayChannelCount2(mDecayBuffer[tPart][0]);
+      tPartBuf->SetDecayChannelCount3(mDecayBuffer[tPart][1]);
+      tPartBuf->SetNumber(tPart);
+      tNum = mDB->AddParticleType(tPartBuf);
+    }
+
+  ifstream in2("i200STAR.m");
+  while (in2.getline(str,200))
+    {
+      tIter=0;
+      //TWO-BODY DECAYS
+      if(str[0] == 's' && str[1] == 'e' && str[2] == '[')
+       {
+         // Reading in the decay channels
+         char *tLBrackert, *tFirstComma, *tSecondComma, *tThirdComma, *tRBracket;
+         char *tFather, *tDaughter1, *tDaughter2, *tBRatio;
+         
+         tLBrackert = strchr(str,'[');
+         tFirstComma = strchr(str,',');
+         tSecondComma = strchr(tFirstComma+1,',');
+         tThirdComma = strchr(tSecondComma+1,',');
+         tRBracket = strchr(tThirdComma,']');
+         if (!(tLBrackert && tFirstComma && tSecondComma && tThirdComma && tRBracket))
+           PRINT_DEBUG_1("Malformed line!: " << str);
+         
+         tFather    = strndup(tLBrackert+1,   tFirstComma-tLBrackert-1);
+         tDaughter1 = strndup(tFirstComma+1,  tSecondComma-tFirstComma-1);
+         tDaughter2 = strndup(tSecondComma+1, tThirdComma-tSecondComma-1);
+         tBRatio    = strndup(tThirdComma+1,  tRBracket-tThirdComma-1);
+         
+         // Getting the ratio
+         char *tMiddle, *tRatComponent;
+         double tRatio = 1.0;
+         
+         tMiddle = strchr(tBRatio,'_');
+         if (!tMiddle) tMiddle = strchr(tBRatio,' ');
+         while(tMiddle)
+           {
+             tRatComponent = strndup(tBRatio, tMiddle-tBRatio);
+             if (strchr(tRatComponent,'/'))
+               tRatio *= atof(tRatComponent)/atof(tRatComponent+2);
+             else
+               tRatio *= atof(tRatComponent);
+             tBRatio = tMiddle+1;              
+             tMiddle = strchr(tBRatio,'_');
+             if (!tMiddle) tMiddle = strchr(tBRatio,' ');
+           }
+         if (strchr(tBRatio,'/'))
+           tRatio *= atof(tBRatio)/atof(tBRatio+2);
+         else
+           tRatio *= atof(tBRatio);
+         
+         DecayChannel *newChannel = new DecayChannel(tRatio, mDB->GetParticleTypeIndex(tDaughter1), mDB->GetParticleTypeIndex(tDaughter2), -1);
+         //      if (mDB->GetParticleType(tDaughter1)->GetMass() + 
+         //          mDB->GetParticleType(tDaughter2)->GetMass()
+         //          < mDB->GetParticleType(tFather)->GetMass()) {
+         //        (mDB->GetParticleType(tFather))->AddDecayChannel(*newChannel);
+         //      }
+         (mDB->GetParticleType(tFather))->AddDecayChannel(*newChannel);
+       }
+      
+      if(str[0] == 's' && str[1] == 'e' && str[2] == '3')
+       {
+         // Reading in the decay channels
+         char *tLBrackert, *tFirstComma, *tSecondComma, *tThirdComma, *tFourthComma, *tRBracket;
+         char *tFather, *tDaughter1, *tDaughter2, *tDaughter3, *tBRatio;
+         
+         tLBrackert = strchr(str,'[');
+         tFirstComma = strchr(str,',');
+         tSecondComma = strchr(tFirstComma+1,',');
+         tThirdComma = strchr(tSecondComma+1,',');
+         tFourthComma = strchr(tThirdComma+1,',');
+         tRBracket = strchr(tThirdComma,']');
+
+         if (!(tLBrackert && tFirstComma && tSecondComma && tThirdComma && tFourthComma && tRBracket))
+           PRINT_DEBUG_1("Malformed line!: " << str);
+         
+         tFather    = strndup(tLBrackert+1,   tFirstComma-tLBrackert-1);
+         tDaughter1 = strndup(tFirstComma+1,  tSecondComma-tFirstComma-1);
+         tDaughter2 = strndup(tSecondComma+1, tThirdComma-tSecondComma-1);
+         tDaughter3 = strndup(tThirdComma+1,  tFourthComma-tThirdComma-1);
+         tBRatio    = strndup(tFourthComma+1, tRBracket-tFourthComma-1);
+         
+         // Getting the ratio
+         char *tMiddle, *tRatComponent;
+         double tRatio = 1.0;
+         
+         tMiddle = strchr(tBRatio,'_');
+         if (!tMiddle) tMiddle = strchr(tBRatio,' ');
+         while(tMiddle)
+           {
+             tRatComponent = strndup(tBRatio, tMiddle-tBRatio);
+             if (strchr(tRatComponent,'/'))
+               tRatio *= atof(tRatComponent)/atof(tRatComponent+2);
+             else
+               tRatio *= atof(tRatComponent);
+             tBRatio = tMiddle+1;              
+             tMiddle = strchr(tBRatio,'_');
+             if (!tMiddle) tMiddle = strchr(tBRatio,' ');
+           }
+         if (strchr(tBRatio,'/'))
+           tRatio *= atof(tBRatio)/atof(tBRatio+2);
+         else
+           tRatio *= atof(tBRatio);
+         
+         DecayChannel *newChannel = new DecayChannel(tRatio, mDB->GetParticleTypeIndex(tDaughter1), mDB->GetParticleTypeIndex(tDaughter2), mDB->GetParticleTypeIndex(tDaughter3));
+         //      if (mDB->GetParticleType(tDaughter1)->GetMass() + 
+         //          mDB->GetParticleType(tDaughter2)->GetMass() +
+         //          mDB->GetParticleType(tDaughter3)->GetMass()
+         //          < mDB->GetParticleType(tFather)->GetMass())
+         (mDB->GetParticleType(tFather))->AddDecayChannel(*newChannel);
+       }
+    }
+
+  ifstream in3("pdgcodes.m");
+  while (in3.getline(str,200))
+    {
+      string tName;
+      int tCode;
+      
+      std::stringstream tIS(str);
+      tIS >> tName >> tCode;
+      mDB->GetParticleType(tName)->SetPDGCode(tCode);
+    }
+  
+}
+
+
+void Parser::ReadShare()
+{
+  char str[50];
+  char str1[200];
+    
+  ParticleType *tPartBuf;
+  int tNum;
+           
+  PRINT_DEBUG_1("Reading from |"<<(mInputDir+"/"+"particles.data").Data()<<"|");
+  ifstream in((mInputDir+"/"+"particles.data").Data());
+    
+  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, I3, I, q, s, aq, as, c, ac, mc;
+         
+         in>>mass>>gamma>>spin>>I>>I3>>q>>s>>aq>>as>>c>>ac>>mc;
+         number++;
+         PRINT_DEBUG_2(number<<" "<<str<<" "<<mass<<" "<<gamma<<" "<<spin<<" "<<I<<" "<<I3<<" "<<q<<" "<<aq<<" "<<s<<" "<<as<<" "<<c<<" "<<ac<<" "<<mc);
+         
+         tPartBuf = new ParticleType();
+         tPartBuf->SetName(str);
+         tPartBuf->SetMass(mass);
+         tPartBuf->SetGamma(gamma);
+         tPartBuf->SetSpin(spin);
+         tPartBuf->SetBarionN((int) ((q+s+c)/3. - (aq+as+ac)/3.) );
+         tPartBuf->SetCharmN((int) (c - ac));
+         tPartBuf->SetStrangeness((int) (as-s));
+         tPartBuf->SetI(I);
+         tPartBuf->SetI3(I3);
+         tPartBuf->SetPDGCode((int) mc);
+         tPartBuf->SetNumber(number);
+         tNum = mDB->AddParticleType(tPartBuf);
+       }
+      in.close();
+    }
+  else 
+    {
+      PRINT_MESSAGE("File "<<(mInputDir+"/"+"particles.data").Data()<<" containing particle data not found!");
+      PRINT_MESSAGE("Please set the correct path to this file in the input parameter file");
+      PRINT_MESSAGE("Aborting!");
+      exit(0);
+    }
+  
+  ifstream in2((mInputDir+"/decays.data").Data());
+  if ((in2) && (in2.is_open()))
+    {
+      //START OF HEAD-LINE
+      in2.ignore(200,'\n');
+      in2.ignore(200,'\n');
+      in2.ignore(200,'\n');
+      //END OF HEAD-LINE
+
+      char tFather[20], tDaughter1[20], tDaughter2[20], tDaughter3[20];
+      double tBRatio, tRatio;
+      int CGcoeff; // complete branching ratio by Clebsch-Gordan coefficient: 0-no 1-yes
+       
+      while (in2>>str)
+       {
+         if (*str == '#')
+           {
+             in2.getline(str1,200);
+             continue;
+           }
+         in2>>tDaughter1>>tDaughter2>>tDaughter3;
+         if (!mDB->ExistsParticleType(tDaughter1)) {
+           PRINT_MESSAGE("Did not find the daughter 1 particle: " << tDaughter1);
+           PRINT_MESSAGE("Not adding channel");
+           in2.getline(str1,200);
+           continue;
+         }
+         if (!mDB->ExistsParticleType(tDaughter2)) {
+           PRINT_MESSAGE("Did not find the daughter 2 particle: " << tDaughter2);
+           PRINT_MESSAGE("Not adding channel");
+           in2.getline(str1,200);
+           continue;
+         }
+         if ((*tDaughter3>65)&&(*tDaughter3<122)&&(!mDB->ExistsParticleType(tDaughter3))) {
+           PRINT_MESSAGE("Did not find the daughter 3 particle: " << tDaughter3);
+           PRINT_MESSAGE("Not adding channel");
+           in2.getline(str1,200);
+           continue;
+         }
+
+         strcpy(tFather,str);
+         PRINT_DEBUG_2(tFather<<"\t"<<tDaughter1<<"\t"<<tDaughter2<<"\t");
+         if ((*tDaughter3>65)&&(*tDaughter3<122)) // check if first char is a letter - if yes then 3-body decay
+           {
+             in2>>tBRatio>>CGcoeff;
+             PRINT_DEBUG_2(tDaughter3<<" (3-body decay)\t");
+             if (mDB->ExistsParticleType(tFather)) {
+               mDB->GetParticleType(tFather)->SetDecayChannelCount3(mDB->GetParticleType(tFather)->GetDecayChannelCount3()+1);
+             
+               tRatio=tBRatio;
+               DecayChannel *newChannel = new DecayChannel(tRatio, mDB->GetParticleTypeIndex(tDaughter1), mDB->GetParticleTypeIndex(tDaughter2), mDB->GetParticleTypeIndex(tDaughter3));
+               if (mDB->GetParticleType(tDaughter1)->GetMass() + 
+                   mDB->GetParticleType(tDaughter2)->GetMass() +
+                   mDB->GetParticleType(tDaughter3)->GetMass()
+                   < mDB->GetParticleType(tFather)->GetMass())
+                 (mDB->GetParticleType(tFather))->AddDecayChannel(*newChannel);
+               
+               tRatio=tBRatio;
+               PRINT_DEBUG_2(tBRatio << '\t' << tRatio);
+             }
+             else {
+               PRINT_MESSAGE("Did not find the father particle: " << tFather);
+               PRINT_MESSAGE("Not adding channel");
+             }
+           }
+         else // 2-body decay      
+           {
+             tBRatio=atof(tDaughter3);
+           
+             in2>>CGcoeff;
+             PRINT_DEBUG_2(" (2-body decay)\t");
+             if (mDB->ExistsParticleType(tFather)) {
+               mDB->GetParticleType(tFather)->SetDecayChannelCount2(mDB->GetParticleType(tFather)->GetDecayChannelCount2()+1);
+               
+               // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+               if (CGcoeff) // complete branching ratio by Clebsch-Gordan coefficient
+                 {
+                   double j1, m1, j2, m2, J, M, CB;
+                   J=mDB->GetParticleType(tFather)->GetI();
+                   M=mDB->GetParticleType(tFather)->GetI3();
+                   j1=mDB->GetParticleType(tDaughter1)->GetI();
+                   m1=mDB->GetParticleType(tDaughter1)->GetI3();
+                   j2=mDB->GetParticleType(tDaughter2)->GetI();
+                   m2=mDB->GetParticleType(tDaughter2)->GetI3();
+                   PRINT_DEBUG_2(" "<<J<<" "<<M<<" "<<j1<<" "<<m1<<" "<<j2<<" "<<m2);
+                   
+                   CB = ClebschGordan(J, M, j1, m1, j2, m2);
+                   tRatio = CB*CB * tBRatio;
+                   
+                   // Multiply the Clebsh by two?
+                   // The same spin, mass, strangeness
+                   // and different I3?
+                   if ((fabs(mDB->GetParticleType(tDaughter1)->GetSpin() - mDB->GetParticleType(tDaughter2)->GetSpin()) < 0.01) &&
+                       (fabs(mDB->GetParticleType(tDaughter1)->GetMass() - mDB->GetParticleType(tDaughter2)->GetMass()) < 0.01) &&
+                       (mDB->GetParticleType(tDaughter1)->GetStrangeness() == mDB->GetParticleType(tDaughter2)->GetStrangeness()) &&
+                       (fabs(mDB->GetParticleType(tDaughter1)->GetI3() - mDB->GetParticleType(tDaughter2)->GetI3()) > 0.01))
+                     {
+                       PRINT_DEBUG_2("Multuplying Clebsch by two for " << tFather << "->" << tDaughter1 << "+" << tDaughter2);
+                       tRatio *= 2.0;
+                     }
+                   
+                   PRINT_DEBUG_2(CB << '\t' << tBRatio << '\t' << tRatio<<"\t"<<CGcoeff);
+                 }
+               
+               else
+                 {
+                   tRatio=tBRatio;
+                   PRINT_DEBUG_2(tBRatio << '\t' << tRatio);
+                 }
+               DecayChannel *newChannel = new DecayChannel(tRatio, mDB->GetParticleTypeIndex(tDaughter1), mDB->GetParticleTypeIndex(tDaughter2), -1);
+               if (mDB->GetParticleType(tDaughter1)->GetMass() + 
+                   mDB->GetParticleType(tDaughter2)->GetMass()
+                   < mDB->GetParticleType(tFather)->GetMass()) 
+                 {
+                   (mDB->GetParticleType(tFather))->AddDecayChannel(*newChannel);
+                   PRINT_DEBUG_2("Added channel " << newChannel << " " << mDB->GetParticleTypeIndex(tFather) << " " << mDB->GetParticleTypeIndex(tDaughter1) << " " << mDB->GetParticleTypeIndex(tDaughter2));
+                 }
+               else 
+                 {
+                   
+                   PRINT_DEBUG_2("Masses do not match! Not adding channel " << newChannel);
+                 }
+             }
+             else {
+               PRINT_MESSAGE("Did not find the father particle: " << tFather);
+               PRINT_MESSAGE("Not adding channel");
+             }
+           }
+       }
+      in2.close();
+    }
+  else {
+    PRINT_MESSAGE("File "<<(mInputDir+"/decays.data").Data()<<" with particle decay channels not found!");
+    PRINT_MESSAGE("No particle decays will be simulated");
+  }
+}
+
+double 
+Parser::ClebschGordan(double aJot, double aEm, double aJot1, double aEm1, double aJot2, double aEm2)
+{
+  int mint, maxt;
+  double cgc = 0.0;
+  int titer;
+  double coef;
+
+  maxt = lrint(aJot1 + aJot2 - aJot);
+  mint = 0;
+  if (lrint(aJot1 - aEm1) < maxt) maxt = lrint(aJot1 - aEm1);
+  if (lrint(aJot2 + aEm2) < maxt) maxt = lrint(aJot2 + aEm2);
+  if (lrint(-(aJot-aJot2+aEm1)) > mint) mint = lrint(-(aJot-aJot2+aEm1));
+  if (lrint(-(aJot-aJot1-aEm2)) > mint) mint = lrint(-(aJot-aJot1-aEm2));
+
+  PRINT_DEBUG_3("mint " << mint << " " <<  aJot1 << " " << aEm1);
+  PRINT_DEBUG_3("maxt " << maxt << " " <<  aJot2 << " " << aEm2);
+
+  for (titer = mint; titer<=maxt; titer ++)
+    {
+      coef = TMath::Power(-1, titer);
+      PRINT_DEBUG_3("coef1 " << coef); 
+      coef *= TMath::Sqrt((2*aJot+1)*
+                         factorials[lrint(aJot1+aEm1)] *
+                         factorials[lrint(aJot1-aEm1)] *
+                         factorials[lrint(aJot2+aEm2)] *
+                         factorials[lrint(aJot2-aEm2)] *
+                         factorials[lrint(aJot+aEm)] *
+                         factorials[lrint(aJot-aEm)]);
+      PRINT_DEBUG_3("coef2 " << coef); 
+      coef /= (factorials[titer] *
+              factorials[lrint(aJot1+aJot2-aJot-titer)] *
+              factorials[lrint(aJot1-aEm1-titer)] *
+              factorials[lrint(aJot2+aEm2-titer)] *
+              factorials[lrint(aJot-aJot2+aEm1+titer)] *
+              factorials[lrint(aJot-aJot1-aEm2+titer)]);
+      PRINT_DEBUG_3("coef3 " << coef); 
+      
+      cgc += coef;
+    }
+
+  cgc *= DeltaJ(aJot1, aJot2, aJot);
+
+  return cgc;
+}
+
+double 
+Parser::DeltaJ(double aJot1, double aJot2, double aJot)
+{
+  double res = TMath::Sqrt(1.0 * 
+                          factorials[lrint(aJot1+aJot2-aJot)] * 
+                          factorials[lrint(aJot1-aJot2+aJot)] * 
+                          factorials[lrint(-aJot1+aJot2+aJot)] / 
+                          factorials[lrint(aJot1+aJot2+aJot+1)]);
+  
+  return res;
+}
+
+void   
+Parser::ReadParameters()
+{
+  // Read the input directory
+  try {
+    mInputDir = sRPInstance->getPar("InputDirSHARE");
+  }
+  catch (STR tError) {
+    PRINT_DEBUG_1("Parser::ReadParameters - Caught exception " << tError);
+    PRINT_MESSAGE("Did not find SHARE input file location.");
+    PRINT_MESSAGE("Using default: '../share'");
+    mInputDir = "../share";
+  }
+}
diff --git a/TTherminator/Therminator/Parser.h b/TTherminator/Therminator/Parser.h
new file mode 100644 (file)
index 0000000..141f7b5
--- /dev/null
@@ -0,0 +1,89 @@
+/******************************************************************************
+ *                      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.                                                                 *
+ *                                                                            *
+ *****************************************************************************/
+#ifndef _BFPW_PARSER_
+#define _BFPW_PARSER_
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+//#include <curses.h>
+#include <iostream>
+#include <fstream>
+#include <map>
+#include <string>
+#include "THGlobal.h"
+#include "ReadPar.h"
+#include "ParticleType.h"
+#include "Integrator.h"
+#include "ParticleDB.h"
+
+class Parser
+{
+ private:
+  char        mNameBuffer[369][20];            //local table with particle names
+  double      mMassBuffer[369];                //local table with particle masses
+  double      mGammaBuffer[369];               //local table with particle gammas      
+  double      mSpinBuffer[369];                //local table with particle spins
+
+  int         mBarionBuffer[369];              //local table with particle barion number
+  int         mStrangeBuffer[369];             //local table with partcle strangeless number
+  double      mI3Buffer[369];          //local table with particle izospin3
+  int         mDecayBuffer[369][2];
+
+  int         mParticleCount;          //number of particles
+  int         mTypeCountBuffer[369];   //number of each particle
+  
+  ParticleDB *mDB;
+
+  TString     mInputDir;
+  
+  double DeltaJ(double aJot1, double aJot2, double aJot);
+  double ClebschGordan(double aJot, double aEm, double aJot1, double aEm1, double aJot2, double aEm2);
+  void   ReadParameters();
+
+ public:
+  Parser(ParticleDB *aDB);
+  ~Parser();
+  void    ReadInput(); //read input1 "tables.m" and "i200STAR.m"
+  void    ReadShare(); //read input from SHaRe format: particles.data & decays.data
+
+  int     GetParticleCount(); 
+       
+  int     check(char *,char *); //check particules by name
+       
+  char   *GetParticleName(int);        
+  double  GetParticleMass(int);
+  double  GetParticleGamma(int);
+  double  GetParticleSpin(int);        
+  int     GetParticleBarionN(int);
+  int     GetParticleStrangeN(int);
+  double  GetParticleI3(int);
+  int     GetParticleDecayChannelCount(int,int);
+  int     GetParticleNumber(int);
+};
+
+#endif
diff --git a/TTherminator/Therminator/Particle.cxx b/TTherminator/Therminator/Particle.cxx
new file mode 100644 (file)
index 0000000..d87cd20
--- /dev/null
@@ -0,0 +1,204 @@
+/******************************************************************************
+ *                      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 "Particle.h"
+#include "ParticleType.h"
+#include <TMath.h>
+#define FIELDWIDTH 15
+
+Particle::Particle()
+{
+}
+
+Particle::Particle(double aRapidity, double aPt, double aPhip, 
+                  double aAlfam, double aRho, double aPhis, double aTau,
+                  ParticleType *aType)
+{
+  mPartType = aType;
+  mDecayed = 0;
+  mHasFather = -1;
+  
+  px = aPt*TMath::Cos(aPhip);
+  py = aPt*TMath::Sin(aPhip);
+  double tMt = TMath::Hypot(GetMass(),aPt);
+  pz = tMt*TMath::SinH(aRapidity);
+  
+  rx = aRho*TMath::Cos(aPhis);
+  ry = aRho*TMath::Sin(aPhis);
+
+// New method of calculating rz, rt
+  rz = aTau*TMath::SinH(aAlfam);
+  rt = aTau*TMath::CosH(aAlfam);
+}
+
+Particle::Particle(ParticleType *aType, 
+                  double aPx, double aPy, double aPz, 
+                  double aRx, double aRy, double aRz,
+                  double aTime)
+{
+  mPartType = aType;
+  mDecayed = 0;
+  mHasFather = -1;
+  
+  px = aPx;
+  py = aPy;
+  pz = aPz;
+
+  rx = aRx;
+  ry = aRy;
+  rz = aRz;
+  rt = aTime;
+}
+
+
+Particle::Particle(const Particle& aParticle)
+{
+  mPartType = aParticle.GetParticleType();
+  mDecayed = 0;
+  mHasFather = aParticle.GetFather();
+  
+  px = aParticle.px;
+  py = aParticle.py;
+  pz = aParticle.pz;
+  
+  rx = aParticle.rx;
+  ry = aParticle.ry;
+  rz = aParticle.rz;
+  rt = aParticle.rt;
+}
+
+Particle::~Particle()
+{
+  /* no-op */
+}
+
+void   
+Particle::WriteParticle(ostream *aOuts)
+{
+  aOuts->flags(ios::left | ios::scientific);
+  aOuts->width(10);
+  (*aOuts) << mPartType->GetPDGCode();
+  aOuts->width(FIELDWIDTH);
+  (*aOuts) << px;
+  aOuts->width(FIELDWIDTH);
+  (*aOuts) << py;
+  aOuts->width(FIELDWIDTH);
+  (*aOuts) << pz;
+  aOuts->width(FIELDWIDTH);
+  (*aOuts) << GetEnergy();
+  aOuts->width(FIELDWIDTH);
+  (*aOuts) << GetMass();
+  aOuts->width(FIELDWIDTH);
+  (*aOuts) << rx;
+  aOuts->width(FIELDWIDTH);
+  (*aOuts) << ry;
+  aOuts->width(FIELDWIDTH);
+  (*aOuts) << rz;
+  aOuts->width(FIELDWIDTH);
+  (*aOuts) << rt;
+//  (*aOuts) << mPartType->GetCharge() << "   ";
+//  if (mHasFather>-1) (*aOuts) << mHasFather << "   ";
+  aOuts->width(6);
+  (*aOuts) << (mHasFather);
+  aOuts->width(2);
+  (*aOuts) << (mDecayed);
+}
+
+
+double 
+Particle::Pt()
+{
+  return TMath::Hypot(px, py);
+}
+
+double 
+Particle::Rapidity()
+{
+  double tE = TMath::Sqrt(px*px+py*py+pz*pz+GetMass()*GetMass());
+  return 0.5*TMath::Log((tE-pz)/(tE+pz));
+}
+
+ParticleType* 
+Particle::GetParticleType() const
+{
+  return mPartType;
+}
+
+int 
+Particle::HadDecayed()
+{
+  return mDecayed;
+}
+
+double 
+Particle::GetMass()
+{
+  return mPartType->GetMass();
+}
+
+double        
+Particle::GetEnergy()
+{
+  return TMath::Sqrt(GetMass()*GetMass()+px*px+py*py+pz*pz);
+}
+
+
+double 
+Particle::GetI3()
+{
+  return mPartType->GetI3();
+}
+
+double        
+Particle::GetBarionN()
+{
+  return mPartType->GetBarionN();
+}
+
+double        
+Particle::GetStrangeness()
+{
+  return mPartType->GetStrangeness();
+}
+
+void
+Particle::SetDecayed()
+{
+  mDecayed = 1;
+}
+
+void          
+Particle::SetFather(int aFather)
+{
+  mHasFather = aFather;
+}
+
+int           
+Particle::GetFather() const
+{
+  return mHasFather;
+}
+
diff --git a/TTherminator/Therminator/Particle.h b/TTherminator/Therminator/Particle.h
new file mode 100644 (file)
index 0000000..74cfca0
--- /dev/null
@@ -0,0 +1,75 @@
+/******************************************************************************
+ *                      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.                                                                 *
+ *                                                                            *
+ *****************************************************************************/
+#ifndef _BFPW_PARTICLE_
+#define _BFPW_PARTICLE_
+
+#include <iostream>
+#include "THGlobal.h"
+
+class ParticleType;
+
+class Particle
+{
+ public:
+  Particle();
+  Particle(double aRapidity, double aPt, double aPhip, 
+          double aAlfam, double aRho, double aPhis, double aTau,
+          ParticleType *aType);
+  Particle(ParticleType *aType, 
+          double aPx, double aPy, double aPz, 
+          double aRx, double aRy, double aRz,
+          double aTime);
+  Particle(const Particle& aParticle);
+  ~Particle();
+
+  double        Pt();
+  double        Rapidity();
+  ParticleType *GetParticleType() const;
+  int           HadDecayed();
+  double        GetMass();
+  double        GetI3();
+  double        GetBarionN();
+  double        GetStrangeness();
+  double        GetEnergy();
+  int           GetFather() const;
+  
+  void          WriteParticle(ostream *aOuts);
+  void          SetDecayed();
+  void          SetFather(int aFather);
+  
+  double px, py, pz;
+  double rx, ry, rz, rt;
+
+ private:
+  int mDecayed;
+  int mHasFather;
+  ParticleType* mPartType;
+};
+
+
+
+#endif
diff --git a/TTherminator/Therminator/ParticleDB.cxx b/TTherminator/Therminator/ParticleDB.cxx
new file mode 100644 (file)
index 0000000..5c2f035
--- /dev/null
@@ -0,0 +1,75 @@
+/******************************************************************************
+ *                      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/0011222,         *
+ * accessibile at: http://www.arxiv.org/nucl-th/0011222                       *
+ *                                                                            *
+ * 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-ex/0011222 and the published version of it, when         *
+ * available.                                                                 *
+ *                                                                            *
+ *****************************************************************************/
+#include "ParticleDB.h"
+
+ParticleDB::ParticleDB()
+{
+  mParticleTable.clear();
+  mParticleNames.clear();
+}
+
+ParticleDB::~ParticleDB()
+{
+}
+
+int           
+ParticleDB::AddParticleType(ParticleType *aPartType)
+{
+  mParticleTable.push_back(*aPartType);
+  mParticleNames[aPartType->GetName()] = mParticleTable.size()-1;
+  return  mParticleTable.size()-1;
+}
+
+ParticleType* 
+ParticleDB::GetParticleType(int aIndex)
+{
+  return &(mParticleTable[aIndex]);
+}
+
+ParticleType* 
+ParticleDB::GetParticleType(std::string aName)
+{
+  return &(mParticleTable[mParticleNames[aName]]);
+}
+
+int           
+ParticleDB::GetParticleTypeIndex(std::string aName)
+{
+  return mParticleNames[aName];
+}
+
+int           
+ParticleDB::GetParticleTypeCount()
+{
+  return mParticleTable.size();
+}
+
+int           
+ParticleDB::ExistsParticleType(std::string aName)
+{
+  return mParticleNames.count(aName);
+}
diff --git a/TTherminator/Therminator/ParticleDB.h b/TTherminator/Therminator/ParticleDB.h
new file mode 100644 (file)
index 0000000..9467ee5
--- /dev/null
@@ -0,0 +1,56 @@
+/******************************************************************************
+ *                      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.                                                                 *
+ *                                                                            *
+ *****************************************************************************/
+#ifndef _BFPW_PARTICLE_DATABASE_
+#define _BFPW_PARTICLE_DATABASE_
+
+#include <map>
+#include <vector>
+#include <string>
+#include "THGlobal.h"
+#include "ParticleType.h"
+#include "DecayChannel.h"
+
+class ParticleDB
+{
+ public:
+  ParticleDB();
+  ~ParticleDB();
+  
+  int           AddParticleType(ParticleType *aPartType);
+  ParticleType* GetParticleType(int aIndex);
+  ParticleType* GetParticleType(std::string aName);
+  int           GetParticleTypeIndex(std::string aName);
+  int           GetParticleTypeCount();
+  int           ExistsParticleType(std::string aName);
+  
+ private:
+  std::vector<ParticleType> mParticleTable;
+  std::map<std::string, int>     mParticleNames;
+};
+
+
+#endif
diff --git a/TTherminator/Therminator/ParticleDecayer.cxx b/TTherminator/Therminator/ParticleDecayer.cxx
new file mode 100644 (file)
index 0000000..a52a9a2
--- /dev/null
@@ -0,0 +1,315 @@
+/******************************************************************************
+ *                      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 "THGlobal.h"
+#include "ParticleDecayer.h"
+#include "DecayChannel.h"
+#include "DecayTable.h"
+#include <TMath.h>
+#include <TDatime.h>
+
+inline Double_t ParticleDecayer::BreitWigner(Double_t Mass, Double_t Gamma) const{
+  Double_t x,y;
+
+  y=mRandom->Rndm();
+  x=Mass+Gamma/2*TMath::Tan(TMath::Pi()*(y-0.5));
+
+  return x;
+}
+
+ParticleDecayer::ParticleDecayer(ParticleDB *aDB)
+{
+  TDatime dat;
+
+  mDB = aDB;
+  mRandom = new TRandom2();
+  mRandom->SetSeed(dat.Get() / 2 * 3);
+}
+
+ParticleDecayer::~ParticleDecayer()
+{
+  delete mRandom;
+}
+
+void
+ParticleDecayer::DecayParticle(Particle *aFather, Particle** aDaughter1, Particle** aDaughter2, Particle** aDaughter3)
+{
+  ParticleType *tType = aFather->GetParticleType();
+#ifdef _RESCALE_CHANNELS_
+  int tChannelIndex = (aFather->GetParticleType())->GetTable()->ChooseDecayChannel(mRandom->Rndm());
+#else
+  Double_t tProb = mRandom->Rndm();
+  int tChannelIndex = (aFather->GetParticleType())->GetTable()->ChooseDecayChannelOrNot(tProb);
+  if (tChannelIndex == -1) {
+    (*aDaughter1) =  NULL;
+    (*aDaughter2) =  NULL;
+    (*aDaughter3) =  NULL; 
+    
+    DecayTable *tab = (aFather->GetParticleType())->GetTable();
+    PRINT_DEBUG_3("Not decaying " << (aFather->GetParticleType())->GetName() << " for prob " << tProb);
+    for (int tIter=0; tIter<=tab->GetChannelCount(); tIter++) {
+      PRINT_DEBUG_3(mDB->GetParticleType(tab->GetDecayChannel(tIter)->GetParticle1())->GetName() << " ");
+      PRINT_DEBUG_3(mDB->GetParticleType(tab->GetDecayChannel(tIter)->GetParticle2())->GetName() << " ");
+      if (tab->GetDecayChannel(tIter)->GetParticle3()>-1)
+       PRINT_DEBUG_3(mDB->GetParticleType(tab->GetDecayChannel(tIter)->GetParticle3())->GetName() << " ");
+      PRINT_DEBUG_3(tab->GetDecayChannel(tIter)->GetBranchingRatio() << endl);
+    }
+    
+    return;
+  }
+
+#endif
+  const DecayChannel *tChan = (aFather->GetParticleType())->GetTable()->GetDecayChannel(tChannelIndex);
+
+  if (tChan->Is3Particle())
+    {
+      // It is a 3-body decay channel
+      ParticleType *tType1 = mDB->GetParticleType(tChan->GetParticle1());
+      ParticleType *tType2 = mDB->GetParticleType(tChan->GetParticle2());
+      ParticleType *tType3 = mDB->GetParticleType(tChan->GetParticle3());
+      
+      double tE = aFather->GetEnergy();
+      double tM = aFather->GetMass();          /*MCH uncommented*/
+      //      double tFatherMass=aFather->GetMass();
+      //      double tFatherGamma=tType->GetGamma();
+      double tM1 = tType1->GetMass();
+      double tM2 = tType2->GetMass();
+      double tM3 = tType3->GetMass();
+/* MCH commented begin
+      double tM;
+      do {
+       tM=BreitWigner(tFatherMass,tFatherGamma);
+      }
+      while (tM1+tM2+tM3>tM);
+MCH commented end*/
+      double tES1, tES2, tP1, tP2, tCos12, tZ;
+      
+      do 
+       {
+         // Generate E1 and E2 with the Momnte-Carlo method
+         do 
+           {
+             tES1 = mRandom->Rndm() * (tM - tM2 - tM3 - tM1) + tM1;
+             tES2 = mRandom->Rndm() * (tM - tM1 - tM3 - tM2) + tM2;
+           }
+         while (tES1+tES2 > tM); // The sum of both energies must be smaller than the resonance mass
+         
+         tP1  = TMath::Sqrt(tES1*tES1 - tM1*tM1);
+         tP2  = TMath::Sqrt(tES2*tES2 - tM2*tM2);
+         
+         tZ = tM - tES1 - tES2;
+         tZ *= tZ;
+         tCos12 = (tZ - tP1*tP1 - tP2*tP2 - tM3*tM3)/(2*tP1*tP2);
+       }
+      while ((tCos12 < -1.0) || (tCos12 > 1.0)); // Cos Theta must exist (be within -1.0 to 1.0 )
+      
+      double tTime;
+      if (tType->GetGamma() == 0.0)
+       tTime = 1.0e10;
+      else {
+       double tTau0 = tE/(tType->GetMass()*tType->GetGamma());
+       // When it decays
+       tTime = -tTau0*TMath::Log(mRandom->Rndm());
+      }
+      
+      // Decay coordinates
+      double rxr = aFather->rx + (aFather->px/tE)*tTime;
+      double ryr = aFather->ry + (aFather->py/tE)*tTime;
+      double rzr = aFather->rz + (aFather->pz/tE)*tTime;
+      double rtr = aFather->rt + tTime;
+      
+      double tPxr2 = tP2 * TMath::Sqrt(1-tCos12*tCos12);
+      double tPzr2 = tP2*tCos12;
+      double tPxr3 = - tPxr2;
+      double tPzr3 = - (tP1 + tPzr2);
+      double tP3 = TMath::Hypot(tPxr3, tPzr3);
+      double tES3 = TMath::Hypot(tM3, tP3);
+
+      // Generating Euler angles
+      double tPhi = mRandom->Rndm() * 2 * TMath::Pi();
+      double tKsi = mRandom->Rndm() * 2 * TMath::Pi();
+      double tCosTh = mRandom->Rndm() * 2.0 - 1.0;
+
+      double sp = TMath::Sin(tPhi);
+      double cp = TMath::Cos(tPhi);
+      double sk = TMath::Sin(tKsi);
+      double ck = TMath::Cos(tKsi);
+      double st = TMath::Sqrt(1.0-tCosTh*tCosTh);
+      double ct = tCosTh;
+
+      // Rotating the whole system
+      double tPxp1 = - st*ck * tP1;
+      double tPyp1 = st*sk * tP1;
+      double tPzp1 = ct * tP1;
+      
+      double tPxp2 = (cp*ct*ck - sp*sk)  * tPxr2 + (-st*ck) * tPzr2;
+      double tPyp2 = (-cp*ct*sk - sp*ck) * tPxr2 + (st*sk)  * tPzr2;
+      double tPzp2 = cp*st               * tPxr2 + ct       * tPzr2;
+
+      double tPxp3 = (cp*ct*ck - sp*sk)  * tPxr3 + (-st*ck) * tPzr3;
+      double tPyp3 = (-cp*ct*sk - sp*ck) * tPxr3 + (st*sk)  * tPzr3;
+      double tPzp3 = cp*st               * tPxr3 + ct       * tPzr3;
+      
+      double tVx = aFather->px/aFather->GetEnergy();
+      double tVy = aFather->py/aFather->GetEnergy();
+      double tVz = aFather->pz/aFather->GetEnergy();
+      
+       tES1 = TMath::Sqrt(tM1*tM1+tPxp1*tPxp1+tPyp1*tPyp1+tPzp1*tPzp1);
+       tES2 = TMath::Sqrt(tM2*tM2+tPxp2*tPxp2+tPyp2*tPyp2+tPzp2*tPzp2);
+       tES3 = TMath::Sqrt(tM3*tM3+tPxp3*tPxp3+tPyp3*tPyp3+tPzp3*tPzp3);
+
+      double tV2 = tVx*tVx + tVy*tVy + tVz*tVz;
+      double tGamma = TMath::Power(1-tV2,-0.5);
+      
+      // Boosting by the parent velocity
+      double tVP = tVx*tPxp1 + tVy*tPyp1 + tVz*tPzp1;
+      double tgvp = (tGamma - 1.0) * (1.0/tV2) * tVP;
+
+      double tPx1 =   tPxp1 + (tgvp + tGamma * tES1) * tVx;
+      double tPy1 =   tPyp1 + (tgvp + tGamma * tES1) * tVy;
+      double tPz1 =   tPzp1 + (tgvp + tGamma * tES1) * tVz;
+  
+      tVP = tVx*tPxp2 + tVy*tPyp2 + tVz*tPzp2;
+      tgvp = (tGamma - 1.0) * (1.0/tV2) * tVP;
+
+      double tPx2 =   tPxp2 + (tgvp + tGamma * tES2) * tVx;
+      double tPy2 =   tPyp2 + (tgvp + tGamma * tES2) * tVy;
+      double tPz2 =   tPzp2 + (tgvp + tGamma * tES2) * tVz;
+  
+      tVP = tVx*tPxp3 + tVy*tPyp3 + tVz*tPzp3;
+      tgvp = (tGamma - 1.0) * (1.0/tV2) * tVP;
+
+      double tPx3 =   tPxp3 + (tgvp + tGamma * tES3) * tVx;
+      double tPy3 =   tPyp3 + (tgvp + tGamma * tES3) * tVy;
+      double tPz3 =   tPzp3 + (tgvp + tGamma * tES3) * tVz;
+      
+      tES1 = TMath::Sqrt(tM1*tM1+tPx1*tPx1+tPy1*tPy1+tPz1*tPz1);
+      tES2 = TMath::Sqrt(tM2*tM2+tPx2*tPx2+tPy2*tPy2+tPz2*tPz2);
+      tES3 = TMath::Sqrt(tM3*tM3+tPx3*tPx3+tPy3*tPy3+tPz3*tPz3);
+      
+      (*aDaughter1) = new Particle(tType1, 
+                                  tPx1, tPy1, tPz1,
+                                  rxr, ryr, rzr, 
+                                  rtr);
+      (*aDaughter2) = new Particle(tType2, 
+                                  tPx2, tPy2, tPz2,
+                                  rxr, ryr, rzr, 
+                                  rtr);
+      (*aDaughter3) = new Particle(tType3, 
+                                  tPx3, tPy3, tPz3,
+                                  rxr, ryr, rzr, 
+                                  rtr);
+
+      aFather->SetDecayed();
+    }
+  else
+    {
+      // It is a regular two-body decay channel
+      ParticleType *tType1 = mDB->GetParticleType(tChan->GetParticle1());
+      ParticleType *tType2 = mDB->GetParticleType(tChan->GetParticle2());
+      
+      double tE = aFather->GetEnergy();
+      double tM = aFather->GetMass();          /*MCH uncommented*/
+      //      double tFatherMass=aFather->GetMass();
+      //      double tFatherGamma=tType->GetGamma();
+      double tM1 = tType1->GetMass();
+      double tM2 = tType2->GetMass();
+/* MCH commented begin
+      double tM;
+      do {
+       tM=BreitWigner(tFatherMass,tFatherGamma);
+      }
+      while (tM1+tM2>tM);
+MCH commented end*/
+
+      double tTime;
+      if (tType->GetGamma() == 0.0)
+       tTime = 1.0e10;
+      else {
+       double tTau0 = tE/(tType->GetMass()*tType->GetGamma());
+       // When it decays
+       tTime = -tTau0*TMath::Log(mRandom->Rndm());
+      }
+      
+      // Decay coordinates
+      double rxr = aFather->rx + (aFather->px/tE)*tTime;
+      double ryr = aFather->ry + (aFather->py/tE)*tTime;
+      double rzr = aFather->rz + (aFather->pz/tE)*tTime;
+      double rtr = aFather->rt + tTime;
+      
+      // Decay energy
+      double tMC1 = (tM*tM - (tM1+tM2)*(tM1+tM2));
+      double tMC2 = (tM*tM - (tM1-tM2)*(tM1-tM2));
+      double tMom = TMath::Sqrt(tMC1*tMC2)/(2*tM);
+      double tPhi = mRandom->Rndm() * 2 * TMath::Pi();
+      double tCosTh = mRandom->Rndm()*2.0-1.0;
+      
+      double tPtr = tMom*TMath::Sqrt(1-tCosTh*tCosTh);
+      double tPxr1 = tPtr*TMath::Cos(tPhi);
+      double tPyr1 = tPtr*TMath::Sin(tPhi);
+      double tPzr1 = tMom*tCosTh;
+      
+      double tVx = aFather->px/aFather->GetEnergy();
+      double tVy = aFather->py/aFather->GetEnergy();
+      double tVz = aFather->pz/aFather->GetEnergy();
+      
+      double tES1 = TMath::Sqrt(tM1*tM1+tPxr1*tPxr1+tPyr1*tPyr1+tPzr1*tPzr1);
+      double tES2 = TMath::Sqrt(tM2*tM2+tPxr1*tPxr1+tPyr1*tPyr1+tPzr1*tPzr1);
+      
+      double tV2 = tVx*tVx + tVy*tVy + tVz*tVz;
+      double tGamma = TMath::Power(1-tV2,-0.5);
+      double tVP = tVx*tPxr1 + tVy*tPyr1 + tVz*tPzr1;
+      double tgvp = (tGamma - 1.0) * (1.0/tV2) * tVP;
+      
+      double tPx1 =   tPxr1 + (tgvp + tGamma * tES1) * tVx;
+      double tPy1 =   tPyr1 + (tgvp + tGamma * tES1) * tVy;
+      double tPz1 =   tPzr1 + (tgvp + tGamma * tES1) * tVz;
+  
+      double tPx2 = - tPxr1 + (-tgvp + tGamma * tES2) * tVx;
+      double tPy2 = - tPyr1 + (-tgvp + tGamma * tES2) * tVy;
+      double tPz2 = - tPzr1 + (-tgvp + tGamma * tES2) * tVz;
+
+      (*aDaughter1) = new Particle(tType1, 
+                                  tPx1, tPy1, tPz1,
+                                  rxr, ryr, rzr, 
+                                  rtr);
+      (*aDaughter2) = new Particle(tType2, 
+                                  tPx2, tPy2, tPz2,
+                                  rxr, ryr, rzr, 
+                                  rtr);
+      (*aDaughter3) = 0;
+
+      aFather->SetDecayed();
+  
+    }
+}
+
+void ParticleDecayer::SeedSet(int aSeed)
+{
+  //  mRandom->SetSeed2(aSeed, aSeed * 11 % 9);
+  mRandom->SetSeed(aSeed);
+}
+
diff --git a/TTherminator/Therminator/ParticleDecayer.h b/TTherminator/Therminator/ParticleDecayer.h
new file mode 100644 (file)
index 0000000..b9feb2f
--- /dev/null
@@ -0,0 +1,51 @@
+/******************************************************************************
+ *                      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.                                                                 *
+ *                                                                            *
+ *****************************************************************************/
+#ifndef _BFPW_PARTICLE_DECAYER_
+#define _BFPW_PARTICLE_DECAYER_
+#include "ParticleDB.h"                                
+#include "ParticleType.h"
+#include "Particle.h"
+#include <TRandom2.h>
+
+class ParticleDecayer
+{
+ public:
+  ParticleDecayer(ParticleDB *aDB);
+  ~ParticleDecayer();
+  
+  void DecayParticle(Particle *aFather, Particle** aDaughter1, Particle** aDaughter2, Particle** aDaughter3);
+  void SeedSet(int aSeed);
+  
+ private:
+  inline Double_t BreitWigner(Double_t Mass, Double_t Gamma) const;
+  ParticleDB *mDB;
+  TRandom2 *mRandom;
+  
+};
+
+
+#endif
diff --git a/TTherminator/Therminator/ParticleType.cxx b/TTherminator/Therminator/ParticleType.cxx
new file mode 100644 (file)
index 0000000..7bd7fa3
--- /dev/null
@@ -0,0 +1,142 @@
+/******************************************************************************
+ *                      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 "ParticleType.h"
+#include "DecayTable.h"
+
+ParticleType::ParticleType()
+{
+  mName = "";
+  mNumber=0;
+  mMass=-1.;
+  mStrangeness=-1;
+  mBarionN=-1;
+  mCharmN=-1;
+  mSpin=-1.;
+  mI=-1.;
+  mI3=-1.;
+  mGamma=-1.;
+  mDecayChannelCount2=0;
+  mDecayChannelCount3=0;
+  mTable = new DecayTable();
+  mPDGCode = 0;
+  mFMax = 0.0;
+}
+
+ParticleType::ParticleType(const ParticleType& aParticleType)
+{
+  mName = aParticleType.GetName();
+  mNumber = aParticleType.GetNumber();
+  mMass = aParticleType.GetMass();
+  mStrangeness = aParticleType.GetStrangeness();
+  mBarionN=aParticleType.GetBarionN();
+  mCharmN=aParticleType.GetCharmN();
+  mSpin=aParticleType.GetSpin();
+  mI=aParticleType.GetI();
+  mI3=aParticleType.GetI3();
+  mGamma=aParticleType.GetGamma();
+  mDecayChannelCount2=aParticleType.GetDecayChannelCount2();
+  mDecayChannelCount3=aParticleType.GetDecayChannelCount3();
+  mPDGCode = aParticleType.GetPDGCode();
+  mFMax = aParticleType.GetFMax();
+  mTable = new DecayTable(*(aParticleType.GetTable()));
+}
+
+ParticleType::~ParticleType()
+{
+  if (mTable)
+    delete mTable;
+}
+
+void ParticleType::WriteParticle(int i)
+{
+    mMass=-1;
+    mStrangeness=-1;
+    mBarionN=-1;
+    mCharmN=-1;
+    mSpin=-1;
+    mName="1";
+}
+
+int    ParticleType::GetNumber() const { return mNumber; }
+float  ParticleType::GetMass() const { return mMass; }
+int    ParticleType::GetStrangeness() const { return mStrangeness; }
+int    ParticleType::GetBarionN() const { return mBarionN; }
+int    ParticleType::GetCharmN() const { return mCharmN; }
+float  ParticleType::GetSpin() const { return mSpin; }
+float  ParticleType::GetI() const { return mI; }
+float  ParticleType::GetI3() const { return mI3; }
+float  ParticleType::GetGamma() const { return mGamma; }
+std::string ParticleType::GetName() const { return mName; }
+int    ParticleType::GetDecayChannelCount2() const { return mDecayChannelCount2; }
+int    ParticleType::GetDecayChannelCount3() const { return mDecayChannelCount3; }
+int    ParticleType::GetCharge() { return int(mI3+(mBarionN+mStrangeness)/2.); }       /*MCH int() added */
+int    ParticleType::GetPDGCode() const { return mPDGCode; }
+
+  
+
+void ParticleType::SetNumber(int aNumber) { mNumber = aNumber; }
+void ParticleType::SetMass(float aMass) { mMass = aMass; }
+void ParticleType::SetStrangeness(int aStrangeness) { mStrangeness = aStrangeness; }
+void ParticleType::SetBarionN(int aBarionN) { mBarionN = aBarionN; }
+void ParticleType::SetCharmN(int aCharmN) { mCharmN = aCharmN; }
+void ParticleType::SetSpin(float aSpin) { mSpin = aSpin; }
+void ParticleType::SetI(float aI) { mI = aI; }
+void ParticleType::SetI3(float aI3) { mI3 = aI3; }
+void ParticleType::SetGamma(float aGamma) { mGamma = aGamma; }
+void ParticleType::SetName(char *aName) { mName = aName; }
+void ParticleType::SetDecayChannelCount2(int aDCCount2) { mDecayChannelCount2 = aDCCount2; }
+void ParticleType::SetDecayChannelCount3(int aDCCount3) { mDecayChannelCount3 = aDCCount3; }
+void ParticleType::SetPDGCode(int aCode) { mPDGCode = aCode; }
+
+DecayTable* 
+ParticleType::GetTable() const
+{
+  if (mTable)
+    return mTable;
+  else
+    return NULL;
+}
+
+void        
+ParticleType::AddDecayChannel(DecayChannel aChannel)
+{
+  if (!mTable)
+    mTable = new DecayTable();
+  mTable->AddDecayChannel(aChannel);
+}
+
+float  
+ParticleType::GetFMax() const
+{
+  return mFMax;
+}
+
+void 
+ParticleType::SetFMax(float aFMax)
+{
+  mFMax = aFMax;
+}
diff --git a/TTherminator/Therminator/ParticleType.h b/TTherminator/Therminator/ParticleType.h
new file mode 100644 (file)
index 0000000..b8f61c1
--- /dev/null
@@ -0,0 +1,101 @@
+/******************************************************************************
+ *                      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.                                                                 *
+ *                                                                            *
+ *****************************************************************************/
+#ifndef _BFPW_PARTICLE_TYPE_
+#define _BFPW_PARTICLE_TYPE_
+
+#include <iostream>
+#include "THGlobal.h"
+
+class DecayChannel;
+class DecayTable;
+
+class ParticleType
+{
+ private:
+
+  int         mNumber;         //particle number
+  float       mMass;           //mass
+  int         mStrangeness;    //strangness
+  int         mBarionN;        //barion number
+  int         mCharmN;         //charm
+  float       mSpin;           //spin
+  float       mI;               //izospin
+  float       mI3;             //izospin3
+  float       mGamma;          //szerokosc polowkowa
+  string      mName;           //particle name
+  int         mPDGCode;
+  int         mDecayChannelCount2;     //number of channels in this case
+  int         mDecayChannelCount3;     //number of channels in this case
+  DecayTable *mTable;
+  float       mFMax;
+
+ public:
+  ParticleType();      //constructor
+  ParticleType(const ParticleType& aParticleType);     //copying constructor
+  ~ParticleType();     //destructor
+
+  int    GetNumber() const;
+  float  GetMass() const;
+  int    GetStrangeness() const;
+  int    GetBarionN() const;
+  int    GetCharmN() const;
+  float  GetSpin() const;
+  float  GetI() const;
+  float  GetI3() const;
+  float  GetGamma() const;
+  string GetName() const;
+  int    GetDecayChannelCount2() const;
+  int    GetDecayChannelCount3() const;
+  int    GetCharge();
+  int    GetPDGCode() const;
+  float  GetFMax() const;
+
+  void SetNumber(int aNumber);
+  void SetMass(float aMass);
+  void SetStrangeness(int aStrangeness);
+  void SetBarionN(int aBarionN);
+  void SetCharmN(int aCharmN);
+  void SetSpin(float aSpin);
+  void SetI(float aI);
+  void SetI3(float aI3);
+  void SetGamma(float aGamma);
+  void SetName(char *aName);
+  void SetDecayChannelCount2(int aDCCount2);
+  void SetDecayChannelCount3(int aDCCount3);
+  void SetPDGCode(int aCode);
+  void SetFMax(float aFMax);
+
+  DecayTable* GetTable() const;
+  void        AddDecayChannel(DecayChannel aChannel);
+  
+  void WriteParticle(int);
+  
+  static int GetParticleTypeNumber(std::string aPartName);
+};
+
+
+#endif
diff --git a/TTherminator/Therminator/ReadPar.cxx b/TTherminator/Therminator/ReadPar.cxx
new file mode 100644 (file)
index 0000000..cb193f2
--- /dev/null
@@ -0,0 +1,115 @@
+/******************************************************************************
+ *                      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 "THGlobal.h"
+#include "ReadPar.h"
+#include <iostream>
+#include <fstream>
+#include <sstream>
+#include <iosfwd>
+
+ReadPar::ReadPar()
+{
+  fname = 0;
+}
+
+ReadPar::ReadPar(const char *aFName)
+{
+  fname = strdup(aFName);
+  readFile(aFName);
+}
+
+int ReadPar::readFile(const char *aFName) throw (int)
+{
+  option read_opt;
+  STR buf_str;
+  char buff[200];
+  char dummy[200];
+
+  std::ifstream       infile(aFName);
+  std::istringstream  *instream;
+  
+  if (!infile.is_open())
+    throw RP_Exception_NoParFile;
+
+  while (!infile.eof())
+    {
+      infile.getline(buff, 200);
+      instream = new std::istringstream(buff);
+      memset(dummy,0,200);
+      *instream >> dummy;
+      
+      PRINT_DEBUG_3("Read " << dummy);;
+      read_opt.keyword = dummy;
+      memset(dummy,0,200);
+      *instream >> dummy;
+
+      PRINT_DEBUG_3("Read " << dummy);
+      if (strstr(dummy,"="))
+       {
+         dummy[0]='\0';
+         
+         memset(dummy,0,200);
+         *instream >> dummy;
+         PRINT_DEBUG_3("Read " << dummy);
+         
+         read_opt.value = dummy;
+         options.push_back(read_opt);
+       }
+    }
+  infile.close();
+
+  return 0; 
+}
+
+int ReadPar::printOptions()
+{
+  VOPT::iterator c;
+
+  for (c=options.begin(); c != options.end(); c++)
+    PRINT_DEBUG_3("Keyword: " << c->keyword << " Value: " << c->value);
+
+  return 0;
+}
+
+STR ReadPar::getPar(const char *name) throw(STR)
+{
+  VOPT::iterator c;
+  STR pname(name);
+
+  for (c=options.begin(); c != options.end(); c++)
+    if (c->keyword == pname)
+      {
+       PRINT_DEBUG_2("Returning value " << c->value << " for keyword " << c->keyword);
+       return c->value;
+      }
+  throw *(new STR(name));
+
+  return TString("");
+}
+
+
+
diff --git a/TTherminator/Therminator/ReadPar.h b/TTherminator/Therminator/ReadPar.h
new file mode 100644 (file)
index 0000000..31a7457
--- /dev/null
@@ -0,0 +1,67 @@
+/******************************************************************************
+ *                      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.                                                                 *
+ *                                                                            *
+ *****************************************************************************/
+#ifndef _THERMINATOR_READPAR_
+#define _THERMINATOR_READPAR_
+#include <iostream>
+#include <string>
+#include <vector>
+#include <exception>
+#include "TString.h"
+
+typedef TString STR;
+
+// Ecxeption values
+#define RP_Exception_UnknownException 0
+#define RP_Exception_NoSuchParamter   1
+#define RP_Exception_NoParFile        2
+
+struct struct_option 
+{
+  STR keyword;
+  STR value;
+};
+
+typedef struct struct_option option;
+typedef std::vector<option> VOPT;
+
+class ReadPar 
+{
+ private:
+  char *fname;
+  VOPT options;
+  
+ public:
+  ReadPar(); // Default constructor
+  ReadPar(const char *aFName);
+  
+  int readFile(const char *aFName) throw(int); 
+  int printOptions();
+  STR getPar(const char *name) throw(STR);
+  
+};
+
+#endif
diff --git a/TTherminator/Therminator/THGlobal.h b/TTherminator/Therminator/THGlobal.h
new file mode 100644 (file)
index 0000000..079ab44
--- /dev/null
@@ -0,0 +1,64 @@
+/******************************************************************************
+ *                      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.                                                                 *
+ *                                                                            *
+ *****************************************************************************/
+#ifndef _CF_GLOBAL_H_
+#define _CF_GLOBAL_H_
+
+// Define global types
+
+using namespace std;
+
+// Define compilation specific variables
+
+#define PRINT_MESSAGE(_mes) cout << _mes << endl;
+
+#if _DEBUG_LEVEL_==0
+#define PRINT_DEBUG_3(_mes) 
+#define PRINT_DEBUG_2(_mes) 
+#define PRINT_DEBUG_1(_mes) 
+#elif _DEBUG_LEVEL_==1
+#define PRINT_DEBUG_3(_mes) 
+#define PRINT_DEBUG_2(_mes) 
+#define PRINT_DEBUG_1(_mes) cerr << _mes << endl;
+#elif _DEBUG_LEVEL_==2
+#define PRINT_DEBUG_3(_mes) 
+#define PRINT_DEBUG_2(_mes) cerr << _mes << endl;
+#define PRINT_DEBUG_1(_mes) cerr << _mes << endl;
+#elif _DEBUG_LEVEL_==3
+#define PRINT_DEBUG_3(_mes) cerr << _mes << endl;
+#define PRINT_DEBUG_2(_mes) cerr << _mes << endl;
+#define PRINT_DEBUG_1(_mes) cerr << _mes << endl;
+#endif
+
+#ifdef _GCC2_
+#define STDIOS ios
+#endif
+
+#ifdef _GCC3_
+#define STDIOS ios_base
+#endif
+
+#endif
diff --git a/TTherminator/libTTherminator.pkg b/TTherminator/libTTherminator.pkg
new file mode 100644 (file)
index 0000000..66ea087
--- /dev/null
@@ -0,0 +1,21 @@
+SRCS= AliGenTherminator.cxx \
+      TTherminator.cxx \
+      Therminator/DecayChannel.cxx \
+      Therminator/DecayTable.cxx \
+      Therminator/Event.cxx \
+      Therminator/Hypersurface.cxx \
+      Therminator/Integrator.cxx \
+      Therminator/Parser.cxx \
+      Therminator/Particle.cxx \
+      Therminator/ParticleDB.cxx \
+      Therminator/ParticleDecayer.cxx \
+      Therminator/ParticleType.cxx \
+      Therminator/ReadPar.cxx
+
+HDRS= $(SRCS:.cxx=.h) Therminator/THGlobal.h
+
+DHDR:=TTherminatorLinkDef.h
+
+EXPORT:=AliGenTherminator.h TTherminator.h
+
+EINCLUDE:= TTherminator/Therminator TTherminator