which included commits to RCS files with non-trunk default branches.
--- /dev/null
+// 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;
+}
--- /dev/null
+#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; }
+
+
+
+
+
+
+
+
--- /dev/null
+# $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)"
+
--- /dev/null
+#! /bin/sh
+
+touch libTTherminator.pkg
+
+make libTTherminator.so
--- /dev/null
+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);
+}
--- /dev/null
+///////////////////////////////////////////////////////////////////////////////
+// //
+// 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;
+
+}
--- /dev/null
+///////////////////////////////////////////////////////////////////////////////
+// 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
--- /dev/null
+#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
--- /dev/null
+/******************************************************************************
+ * 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;
+}
--- /dev/null
+/******************************************************************************
+ * 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
--- /dev/null
+/******************************************************************************
+ * 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;
+}
+
--- /dev/null
+/******************************************************************************
+ * 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
--- /dev/null
+/******************************************************************************
+ * 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;
+ }
+
+}
+
--- /dev/null
+/******************************************************************************
+ * 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
--- /dev/null
+#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);
+}
--- /dev/null
+/***************************************************************************
+ * 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_ */
--- /dev/null
+/******************************************************************************
+ * 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);
+}
+
--- /dev/null
+/******************************************************************************
+ * 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
+
--- /dev/null
+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 $^
--- /dev/null
+/******************************************************************************
+ * 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";
+ }
+}
--- /dev/null
+/******************************************************************************
+ * 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
--- /dev/null
+/******************************************************************************
+ * 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;
+}
+
--- /dev/null
+/******************************************************************************
+ * 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
--- /dev/null
+/******************************************************************************
+ * 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);
+}
--- /dev/null
+/******************************************************************************
+ * 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
--- /dev/null
+/******************************************************************************
+ * 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);
+}
+
--- /dev/null
+/******************************************************************************
+ * 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
--- /dev/null
+/******************************************************************************
+ * 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;
+}
--- /dev/null
+/******************************************************************************
+ * 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
--- /dev/null
+/******************************************************************************
+ * 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("");
+}
+
+
+
--- /dev/null
+/******************************************************************************
+ * 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
--- /dev/null
+/******************************************************************************
+ * 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
--- /dev/null
+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