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

/////////////////////////////////////////////////////////////////////////////
// TUHKM is an interface to 
// UHKM 3.0                                                                //
// ( only the HYDJET++ part is implemented)                                // 
// temporary link:                                                         //
// http://lav01.sinp.msu.ru/~igor/hydjet++/hydjet++.txt                    //  
// The main UHKM options are accessable  through this interface.           //
// Uses the TUHKMgen implementation of TGenerator.                         //
// Author of the first implementation: Sergey Zaporozhets                  //
// (zaporozh@sunhe.jinr.ru)                                                // 
// Futhers modifications were made by                                      //
// Ionut Cristian Arsene (i.c.arsene@fys.uio.no)                           //
// & Malinina Liudmila(malinina@lav01.sinp.msu.ru) using as an example     //
//  AliGenTherminator.cxx created by Adam Kisiel                           //   
//                                                                         //
////////////////////////////////////////////////////////////////////////////

#ifndef TUHKMGEN_H
#include "TUHKMgen.h"
#endif
#include "TObjArray.h"
#include "TParticle.h"
#include "Particle.h"
#include "TROOT.h"
#include "TError.h"
#include "TClonesArray.h"
#include "TSystem.h"

#include <string>
#include <iostream>
using namespace std;

//class TObjArray;

ClassImp(TUHKMgen)

TUHKMgen::TUHKMgen() : 
  TGenerator("UHKM","UHKM"),
  fInitialState(0x0),
  fNPprim(0),
  fNPsec(0),
  fHydjetParams(),
  fStableFlagged(0)
{
  // default constructor setting reasonable defaults for initial parameters (central Pb+Pb at 5.5 TeV)

  ParticleAllocator fAllocator;
  List_t fSecondariesList;

  // Set reasonable default values for LHC
  
  fHydjetParams.fSqrtS=5500; //LHC
  fHydjetParams.fAw=207;//Au-Au
  fHydjetParams.fIfb    = 1;      // flag of type of centrality generation (=0 is fixed by fBfix, not 0 
                                   //impact parameter is generated in each event between fBfmin 
                                   //and fBmax according with Glauber model (f-la 30)
  fHydjetParams.fBfix=0.;
  fHydjetParams.fBmin=0.;
  fHydjetParams.fBmax=0.5; //0-5% centrality
  fHydjetParams.fT = 0.170;
  fHydjetParams.fMuB = 0.0;
  fHydjetParams.fMuS = 0.0;
  fHydjetParams.fMuI3 = 0.0;
  fHydjetParams.fThFO = 0.130;
  fHydjetParams.fMu_th_pip = 0.0;
  fHydjetParams.fSeed=0;
  fHydjetParams.fTau=10.;
  fHydjetParams.fSigmaTau=3.;
  fHydjetParams.fR=11.;
  fHydjetParams.fYlmax=4.0;
  fHydjetParams.fUmax=1.1;
  fHydjetParams.fDelta=0.;
  fHydjetParams.fEpsilon=0.;
  fHydjetParams.fWeakDecay=0; //>=0 on ,-1 off
  fHydjetParams.fEtaType=1;//gaus
  fHydjetParams.fCorrS=1.;
  fHydjetParams.fNhsel=2;
  fHydjetParams.fIshad=1;
  fHydjetParams.fPtmin=7.0;
  fHydjetParams.fT0=0.8;
  fHydjetParams.fTau0=0.1;
  fHydjetParams.fNf=0;
  fHydjetParams.fIenglu=0;
  fHydjetParams.fIanglu=0;


  // Set reasonable default values for RHIC
/*
  fHydjetParams.fSqrtS=200; //RHIC          
  fHydjetParams.fAw=197;//Au-Au        
  fHydjetParams.fIfb    = 1;      
  fHydjetParams.fBfix=0.;       
  fHydjetParams.fBmin=0.;       
  fHydjetParams.fBmax=0.5; //0-5% centrality
  fHydjetParams.fT = 0.165;
  fHydjetParams.fMuB = 0.0285;
  fHydjetParams.fMuS = 0.007;
  fHydjetParams.fMuI3 = -0.001;
  fHydjetParams.fThFO = 0.100;
  fHydjetParams.fMu_th_pip = 0.053;
  fHydjetParams.fSeed=0;  
  fHydjetParams.fTau=8.;
  fHydjetParams.fSigmaTau=2.;
  fHydjetParams.fR=10.;
  fHydjetParams.fYlmax=3.3;
  fHydjetParams.fUmax=1.1;
  fHydjetParams.fDelta=0.;
  fHydjetParams.fEpsilon=0.;
  fHydjetParams.fWeakDecay=0;//>=0 on ,-1 off 
  fHydjetParams.fEtaType=1;//gaus
  fHydjetParams.fCorrS=1.;  
  fHydjetParams.fNhsel=2;
  fHydjetParams.fIshad=1;
  fHydjetParams.fPtmin=3.4;
  fHydjetParams.fT0=0.3;
  fHydjetParams.fTau0=0.4;
  fHydjetParams.fNf=2;      
  fHydjetParams.fIenglu=0;     
  fHydjetParams.fIanglu=0;  
*/
   
  strcpy(fParticleFilename, Form("%s/TUHKMgen/UHKM/particles.data", gSystem->Getenv("ALICE_ROOT")));
  strcpy(fDecayFilename, Form("%s/TUHKMgen/UHKM/tabledecay.txt", gSystem->Getenv("ALICE_ROOT")));
  for(Int_t i=0; i<500; i++) {
    fStableFlagPDG[i] = 0;
    fStableFlagStatus[i] = kFALSE;
  }
  fStableFlagged = 0;
  //  cout << "TUHKMgen::TUHKMgen() OUT" << endl;
}

//______________________________________________________________________________
TUHKMgen::~TUHKMgen()
{
  // destructor, deletes the InitialStateHydjet object
  if(fInitialState)
    delete fInitialState;
}

void TUHKMgen::SetAllParametersRHIC()
{
  // Set reasonable input parameters for 0-5% central Au+Au collisions
  // at 200 GeV at RHIC
  SetEcms(200.0);                // RHIC top energy
  SetAw(197);                    // Au+Au
  SetIfb(1);
  SetBfix(0.);
  SetBmin(0.0);                  // 0%
  SetBmax(0.5);                  // 5%
  SetChFrzTemperature(0.165);    // T_ch = 165 MeV
  SetMuB(0.0285);                // mu_B = 28.5 MeV
  SetMuS(0.007);                 // mu_S = 7 MeV
  SetMuQ(-0.001);                // mu_Q = -1 MeV
  SetThFrzTemperature(0.100);    // T_th = 100 MeV
  SetMuPionThermal(0.053);       // mu_th_pion = 53 MeV
  SetSeed(0);                    // use UNIX time
  SetTauB(8.0);                  // tau = 8 fm/c
  SetSigmaTau(2.0);              // sigma_tau = 2 fm/c
  SetRmaxB(10.0);                // fR = 10 fm
  SetYlMax(3.3);                 // fYmax = 3.3
  SetEtaRMax(1.1);               // Umax = 1.1
  SetMomAsymmPar(0.0);           // delta = 0.0
  SetCoordAsymmPar(0.0);         // epsilon = 0.0
  SetFlagWeakDecay(0.0);           // weak decay on (<0 off !!!)
  SetEtaType(1);                 // gaus distributed with fYmax dispersion (0 means boost invariant)
  SetGammaS(1.0);                // gammaS = 1.0 (no strangeness canonical suppresion)
  SetPyquenNhsel(2);             // hydro on, jets on, jet quenching on
  SetPyquenShad(1);              // shadowing on (0 off)
  SetPyquenPtmin(3.4);           // ptmin = 3.4 GeV/c
  SetPyquenT0(0.3);              // T0 = 300 MeV
  SetPyquenTau0(0.4);            // tau0 = 0.4 fm/c
  SetPyquenNf(2);                // 2 flavours
  SetPyquenIenglu(0);            // radiative and collisional energy loss
  SetPyquenIanglu(0);            // small gluon angular distribution
}

void TUHKMgen::SetAllParametersLHC()
{
  // Set reasonable input parameters for 0-5% Pb+Pb collisions at 5.5 TeV at LHC
  SetEcms(5500.0);               // LHC
  SetAw(207);                    // Pb+Pb
  SetIfb(1);
  SetBfix(0.);                  // 0
  SetBmin(0.0);                  // 0%
  SetBmax(0.5);                  // 5%
  SetChFrzTemperature(0.170);    // T_ch = 170 MeV
  SetMuB(0.0);                   // mu_B = 0 MeV
  SetMuS(0.0);                   // mu_S = 0 MeV
  SetMuQ(0.0);                   // mu_Q = 0 MeV
  SetThFrzTemperature(0.130);    // T_th = 130 MeV
  SetMuPionThermal(0.0);         // mu_th_pion = 0 MeV
  SetSeed(0);                    // use UNIX time
  SetTauB(10.0);                 // tau = 10 fm/c
  SetSigmaTau(3.0);              // sigma_tau = 3 fm/c
  SetRmaxB(11.0);                // fR = 11 fm
  SetYlMax(4.0);                 // fYmax = 4.0
  SetEtaRMax(1.1);               // Umax = 1.1
  SetMomAsymmPar(0.0);           // delta = 0.0
  SetCoordAsymmPar(0.0);         // epsilon = 0.0
  SetFlagWeakDecay(0);           // weak decay on (<0 off !!!)
  SetEtaType(1);                 // gaus distributed with fYmax dispersion (0 means boost invariant)
  SetGammaS(1.0);                // gammaS = 1.0 (no strangeness canonical suppresion)
  SetPyquenNhsel(2);             // hydro on, jets on, jet quenching on
  SetPyquenShad(1);              // shadowing on (0 off)
  SetPyquenPtmin(7.0);           // ptmin = 7.0 GeV/c
  SetPyquenT0(0.8);              // T0 = 800 MeV
  SetPyquenTau0(0.1);            // tau0 = 0.4 fm/c
  SetPyquenNf(0);                // 0 flavours
  SetPyquenIenglu(0);            // radiative and collisional energy loss
  SetPyquenIanglu(0);            // small gluon angular distribution
}

TObjArray* TUHKMgen::ImportParticles(const Option_t *)
{
  // Function overloading the TGenerator::ImportParticles() member function.
  // The particles from the local particle list (fSecondariesList) are
  // forwarded to the TGenerator::fParticles 
  fParticles->Clear();
  Int_t nump = 0;
  LPIT_t it,e;
   
  for(it = fSecondariesList.begin(), e = fSecondariesList.end(); it != e; it++) {
    TVector3 pos(it->Pos().Vect());
    TVector3 mom(it->Mom().Vect());
    Float_t m1 = it->TableMass();
    Int_t im1 = it->GetMother();
    Int_t im2 = -1;
    Int_t id1 = -1;
    Int_t id2 = -1;

    Int_t type = it->GetType();  // 0-hydro, 1-jets

    if (im1> -1) {
      TParticle *mother = (TParticle*) (fParticles->UncheckedAt(im1+1));	   
      mother->SetLastDaughter(nump);
      if (mother->GetFirstDaughter()==-1)
 	mother->SetFirstDaughter(nump+1);
    }

    nump++;

    TParticle* p = new TParticle(it->Encoding(), type,	                        //pdg,stat
				 im1, im2, id1, id2,				                //m1,m2,d1,d2
				 mom[0], mom[1], mom[2], TMath::Sqrt(mom.Mag2() + m1 * m1),	//px,py,pz,e
				 pos[0]*1.e-13, pos[1]*1.e-13, pos[2]*1.e-13, 
				 it->T()*1.e-13/3e10);    		        //x,y,z,t
  
     p->SetUniqueID(nump);
     fParticles->Add(p);
  } //end for
 
  fAllocator.FreeList(fSecondariesList);

  return fParticles;
}

Int_t TUHKMgen::ImportParticles(TClonesArray *particles, const Option_t* option)
{
  // Function overloading the TGenerator::ImportParticles() member function.
  // The particles from the local particle list (fSecondariesList) are
  // forwarded to the TGenerator::fParticles

  
  if(particles==0) return 0;
  TClonesArray &particlesR=*particles;
  particlesR.Clear();
 
  Int_t numprim,numsec;  numprim=numsec=0;
  Int_t nump = 0;
  LPIT_t it,e;
  
  for(it = fSecondariesList.begin(), e = fSecondariesList.end(); it != e; it++) {
    TVector3 pos(it->Pos().Vect());  
    TVector3 mom(it->Mom().Vect());  
    Float_t m1 = it->TableMass();  
    Int_t im1 = it->GetMother();  
    Int_t im2 = -1;
    Int_t id1 = -1;
    Int_t id2 = -1;
    
    Int_t type = it->GetType();  
      
    if (im1> -1) {
     // particle not a primary -> set the daughter indexes for the mother particle"<< endl;
      TParticle *mother = (TParticle*) (particlesR.UncheckedAt(im1));
      mother->SetLastDaughter(nump);
      if(mother->GetFirstDaughter()==-1)
	mother->SetFirstDaughter(nump);
    }
      
    
    new (particlesR[nump]) TParticle(it->Encoding(), type,			                        //pdg,stat
				     im1, im2, id1, id2,				                //m1,m2,d1,d2
				     mom[0], mom[1], mom[2], TMath::Sqrt(mom.Mag2() + m1 * m1),	//px,py,pz,e
				     pos[0]*1.e-13, pos[1]*1.e-13, pos[2]*1.e-13, 
				     it->T()*1.e-13/3e10);    		        //x,y,z,t
    
    particlesR[nump]->SetUniqueID(nump);
    nump++;
    numsec++;
  }//end for
  
  fSecondariesList.clear();
  printf("Scan and add prim %d sec %d and all %d particles\n",
  	 numprim,numsec,nump);
  return nump;
}

//______________________________________________________________________________
void TUHKMgen::Initialize()
{
  // Function overloading the TGenerator::Initialize() member function.
  // The Monte-Carlo model is initialized (input parameters are transmitted, 
  // particle list and decay channels are loaded, average multiplicities are calculated, etc.)
 
  fInitialState = new InitialStateHydjet();  
  SetAllParameters();  
  fInitialState->LoadPDGInfo();  
  // set the stable flags
  for(Int_t i=0; i<fStableFlagged; i++) 
    fInitialState->SetPDGParticleStable(fStableFlagPDG[i], fStableFlagStatus[i]); 

  if(!fInitialState->MultIni()) 
    Error("TUHKMgen::Initialize()", "Bad status return from MultIni(). Check it out!! \n");  //
}

void TUHKMgen::Print(const Option_t*) const
{
  cout << "TUHKMgen::Print() method not implemented yet!!" << endl;
}

void TUHKMgen::GenerateEvent()
{
  // Member function overloading the TGenerator::GenerateEvent()
  // The HYDJET++ model is run and the particle lists (fSourceList and fSecondariesList) are filled
  
  fInitialState->Initialize(fSecondariesList, fAllocator);  
 
  if(fSecondariesList.empty())Error("TUHKMgen::GenerateEvent()", "Source particle list empty after fireball initialization!! \n");  //

  // Run the decays

  if(fInitialState->RunDecays())  
    fInitialState->Evolve(fSecondariesList, fAllocator, fInitialState->GetWeakDecayLimit());  
  
}

void TUHKMgen::SetAllParameters() {
  // forward all model input parameters to the InitialStateHydjet object
  // which will handle all the Monte-Carlo simulation using HYDJET++ model

  fInitialState->fParams.fSqrtS = fHydjetParams.fSqrtS;
  fInitialState->fParams.fAw = fHydjetParams.fAw;
  fInitialState->fParams.fIfb = fHydjetParams.fIfb ;
  fInitialState->fParams.fBfix = fHydjetParams.fBfix;
  fInitialState->fParams.fBmin = fHydjetParams.fBmin;
  fInitialState->fParams.fBmax = fHydjetParams.fBmax;
  fInitialState->fParams.fSeed = fHydjetParams.fSeed;
  fInitialState->fParams.fT = fHydjetParams.fT;
  fInitialState->fParams.fMuB = fHydjetParams.fMuB;
  fInitialState->fParams.fMuS = fHydjetParams.fMuS;
  fInitialState->fParams.fMuI3 = fHydjetParams.fMuI3;
  fInitialState->fParams.fThFO = fHydjetParams.fThFO;
  fInitialState->fParams.fMu_th_pip = fHydjetParams.fMu_th_pip;

  fInitialState->fParams.fTau = fHydjetParams.fTau;
  fInitialState->fParams.fSigmaTau = fHydjetParams.fSigmaTau;
  fInitialState->fParams.fR = fHydjetParams.fR;
  fInitialState->fParams.fYlmax = fHydjetParams.fYlmax;
  fInitialState->fParams.fUmax = fHydjetParams.fUmax;
  fInitialState->fParams.fDelta = fHydjetParams.fDelta;
  fInitialState->fParams.fEpsilon = fHydjetParams.fEpsilon;

  fInitialState->fParams.fWeakDecay = fHydjetParams.fWeakDecay;
  fInitialState->fParams.fDecay = 1;                   // run decays
  fInitialState->fParams.fCorrS = fHydjetParams.fCorrS;
  fInitialState->fParams.fEtaType = fHydjetParams.fEtaType;
 
  fInitialState->fParams.fPtmin = fHydjetParams.fPtmin;
  fInitialState->fParams.fNhsel = fHydjetParams.fNhsel;
  fInitialState->fParams.fIshad = fHydjetParams.fIshad;
  fInitialState->fParams.fT0 = fHydjetParams.fT0;
  fInitialState->fParams.fTau0 = fHydjetParams.fTau0;
  fInitialState->fParams.fNf = fHydjetParams.fNf;      
  fInitialState->fParams.fIenglu = fHydjetParams.fIenglu;     
  fInitialState->fParams.fIanglu = fHydjetParams.fIanglu;  

  fInitialState->SetPDGParticleFilename(fParticleFilename);
  fInitialState->SetPDGDecayFilename(fDecayFilename);
  //  fInitialState->SetUseCharmParticles(fUseCharmParticles);
  //  fInitialState->SetWidthRange(fMinWidth, fMaxWidth);
  //  fInitialState->SetMassRange(fMinMass, fMaxMass);
  //  for(Int_t i=0; i<fStableFlagged; i++) 
  //    fInitialState->SetPDGParticleStable(fStableFlagPDG[i], fStableFlagStatus[i]);
  //  cout << "TUHKMgen::SetAllParameters() OUT" << endl;
}
Commit	Line	Data
b1c2e580	1	/////////////////////////////////////////////////////////////////////////////
	2	// TUHKM is an interface to
	3	// UHKM 3.0 //
	4	// ( only the HYDJET++ part is implemented) //
	5	// temporary link: //
	6	// http://lav01.sinp.msu.ru/~igor/hydjet++/hydjet++.txt //
	7	// The main UHKM options are accessable through this interface. //
	8	// Uses the TUHKMgen implementation of TGenerator. //
	9	// Author of the first implementation: Sergey Zaporozhets //
	10	// (zaporozh@sunhe.jinr.ru) //
	11	// Futhers modifications were made by //
	12	// Ionut Cristian Arsene (i.c.arsene@fys.uio.no) //
	13	// & Malinina Liudmila(malinina@lav01.sinp.msu.ru) using as an example //
	14	// AliGenTherminator.cxx created by Adam Kisiel //
	15	// //
	16	////////////////////////////////////////////////////////////////////////////
	17
7b7936e9	18	#ifndef TUHKMGEN_H
b1c2e580	19	#include "TUHKMgen.h"
	20	#endif
	21	#include "TObjArray.h"
	22	#include "TParticle.h"
	23	#include "Particle.h"
	24	#include "TROOT.h"
	25	#include "TError.h"
	26	#include "TClonesArray.h"
	27	#include "TSystem.h"
	28
	29	#include <string>
	30	#include <iostream>
	31	using namespace std;
	32
7b7936e9	33	//class TObjArray;
7b7936e9	34
b1c2e580	35	ClassImp(TUHKMgen)
	36
	37	TUHKMgen::TUHKMgen() :
	38	TGenerator("UHKM","UHKM"),
	39	fInitialState(0x0),
786056a2	40	fNPprim(0),
	41	fNPsec(0),
	42	fHydjetParams(),
7b7936e9	43	fStableFlagged(0)
b1c2e580	44	{
7b7936e9	45	// default constructor setting reasonable defaults for initial parameters (central Pb+Pb at 5.5 TeV)
7b7936e9	46
3fa37a65	47	ParticleAllocator fAllocator;
3fa37a65	48	List_t fSecondariesList;
b1c2e580	49
	50	// Set reasonable default values for LHC
	51
	52	fHydjetParams.fSqrtS=5500; //LHC
	53	fHydjetParams.fAw=207;//Au-Au
	54	fHydjetParams.fIfb = 1; // flag of type of centrality generation (=0 is fixed by fBfix, not 0
	55	//impact parameter is generated in each event between fBfmin
	56	//and fBmax according with Glauber model (f-la 30)
	57	fHydjetParams.fBfix=0.;
	58	fHydjetParams.fBmin=0.;
	59	fHydjetParams.fBmax=0.5; //0-5% centrality
	60	fHydjetParams.fT = 0.170;
	61	fHydjetParams.fMuB = 0.0;
	62	fHydjetParams.fMuS = 0.0;
	63	fHydjetParams.fMuI3 = 0.0;
	64	fHydjetParams.fThFO = 0.130;
	65	fHydjetParams.fMu_th_pip = 0.0;
	66	fHydjetParams.fSeed=0;
	67	fHydjetParams.fTau=10.;
	68	fHydjetParams.fSigmaTau=3.;
	69	fHydjetParams.fR=11.;
	70	fHydjetParams.fYlmax=4.0;
	71	fHydjetParams.fUmax=1.1;
	72	fHydjetParams.fDelta=0.;
	73	fHydjetParams.fEpsilon=0.;
	74	fHydjetParams.fWeakDecay=0; //>=0 on ,-1 off
	75	fHydjetParams.fEtaType=1;//gaus
	76	fHydjetParams.fCorrS=1.;
	77	fHydjetParams.fNhsel=2;
	78	fHydjetParams.fIshad=1;
	79	fHydjetParams.fPtmin=7.0;
	80	fHydjetParams.fT0=0.8;
	81	fHydjetParams.fTau0=0.1;
	82	fHydjetParams.fNf=0;
	83	fHydjetParams.fIenglu=0;
	84	fHydjetParams.fIanglu=0;
	85
	86
	87	// Set reasonable default values for RHIC
	88	/*
	89	fHydjetParams.fSqrtS=200; //RHIC
	90	fHydjetParams.fAw=197;//Au-Au
	91	fHydjetParams.fIfb = 1;
	92	fHydjetParams.fBfix=0.;
	93	fHydjetParams.fBmin=0.;
	94	fHydjetParams.fBmax=0.5; //0-5% centrality
	95	fHydjetParams.fT = 0.165;
	96	fHydjetParams.fMuB = 0.0285;
	97	fHydjetParams.fMuS = 0.007;
	98	fHydjetParams.fMuI3 = -0.001;
	99	fHydjetParams.fThFO = 0.100;
	100	fHydjetParams.fMu_th_pip = 0.053;
	101	fHydjetParams.fSeed=0;
	102	fHydjetParams.fTau=8.;
	103	fHydjetParams.fSigmaTau=2.;
	104	fHydjetParams.fR=10.;
	105	fHydjetParams.fYlmax=3.3;
	106	fHydjetParams.fUmax=1.1;
	107	fHydjetParams.fDelta=0.;
	108	fHydjetParams.fEpsilon=0.;
	109	fHydjetParams.fWeakDecay=0;//>=0 on ,-1 off
	110	fHydjetParams.fEtaType=1;//gaus
	111	fHydjetParams.fCorrS=1.;
	112	fHydjetParams.fNhsel=2;
113	fHydjetParams.fIshad=1;
114	fHydjetParams.fPtmin=3.4;
115	fHydjetParams.fT0=0.3;
116	fHydjetParams.fTau0=0.4;
117	fHydjetParams.fNf=2;
118	fHydjetParams.fIenglu=0;
119	fHydjetParams.fIanglu=0;
120	*/
121
122	strcpy(fParticleFilename, Form("%s/TUHKMgen/UHKM/particles.data", gSystem->Getenv("ALICE_ROOT")));
123	strcpy(fDecayFilename, Form("%s/TUHKMgen/UHKM/tabledecay.txt", gSystem->Getenv("ALICE_ROOT")));
124	for(Int_t i=0; i<500; i++) {
125	fStableFlagPDG[i] = 0;
126	fStableFlagStatus[i] = kFALSE;
127	}
128	fStableFlagged = 0;
7b7936e9	129	// cout << "TUHKMgen::TUHKMgen() OUT" << endl;
b1c2e580	130	}
	131
	132	//______________________________________________________________________________
	133	TUHKMgen::~TUHKMgen()
	134	{
7b7936e9	135	// destructor, deletes the InitialStateHydjet object
b1c2e580	136	if(fInitialState)
b1c2e580	137	delete fInitialState;
b1c2e580	138	}
	139
	140	void TUHKMgen::SetAllParametersRHIC()
	141	{
7b7936e9	142	// Set reasonable input parameters for 0-5% central Au+Au collisions
7b7936e9	143	// at 200 GeV at RHIC
b1c2e580	144	SetEcms(200.0); // RHIC top energy
	145	SetAw(197); // Au+Au
	146	SetIfb(1);
	147	SetBfix(0.);
	148	SetBmin(0.0); // 0%
	149	SetBmax(0.5); // 5%
	150	SetChFrzTemperature(0.165); // T_ch = 165 MeV
	151	SetMuB(0.0285); // mu_B = 28.5 MeV
	152	SetMuS(0.007); // mu_S = 7 MeV
	153	SetMuQ(-0.001); // mu_Q = -1 MeV
	154	SetThFrzTemperature(0.100); // T_th = 100 MeV
	155	SetMuPionThermal(0.053); // mu_th_pion = 53 MeV
	156	SetSeed(0); // use UNIX time
	157	SetTauB(8.0); // tau = 8 fm/c
	158	SetSigmaTau(2.0); // sigma_tau = 2 fm/c
	159	SetRmaxB(10.0); // fR = 10 fm
	160	SetYlMax(3.3); // fYmax = 3.3
	161	SetEtaRMax(1.1); // Umax = 1.1
	162	SetMomAsymmPar(0.0); // delta = 0.0
	163	SetCoordAsymmPar(0.0); // epsilon = 0.0
	164	SetFlagWeakDecay(0.0); // weak decay on (<0 off !!!)
	165	SetEtaType(1); // gaus distributed with fYmax dispersion (0 means boost invariant)
	166	SetGammaS(1.0); // gammaS = 1.0 (no strangeness canonical suppresion)
	167	SetPyquenNhsel(2); // hydro on, jets on, jet quenching on
	168	SetPyquenShad(1); // shadowing on (0 off)
	169	SetPyquenPtmin(3.4); // ptmin = 3.4 GeV/c
	170	SetPyquenT0(0.3); // T0 = 300 MeV
	171	SetPyquenTau0(0.4); // tau0 = 0.4 fm/c
	172	SetPyquenNf(2); // 2 flavours
	173	SetPyquenIenglu(0); // radiative and collisional energy loss
	174	SetPyquenIanglu(0); // small gluon angular distribution
	175	}
	176
	177	void TUHKMgen::SetAllParametersLHC()
	178	{
7b7936e9	179	// Set reasonable input parameters for 0-5% Pb+Pb collisions at 5.5 TeV at LHC
b1c2e580	180	SetEcms(5500.0); // LHC
	181	SetAw(207); // Pb+Pb
	182	SetIfb(1);
	183	SetBfix(0.); // 0
	184	SetBmin(0.0); // 0%
	185	SetBmax(0.5); // 5%
	186	SetChFrzTemperature(0.170); // T_ch = 170 MeV
	187	SetMuB(0.0); // mu_B = 0 MeV
	188	SetMuS(0.0); // mu_S = 0 MeV
	189	SetMuQ(0.0); // mu_Q = 0 MeV
	190	SetThFrzTemperature(0.130); // T_th = 130 MeV
	191	SetMuPionThermal(0.0); // mu_th_pion = 0 MeV
	192	SetSeed(0); // use UNIX time
	193	SetTauB(10.0); // tau = 10 fm/c
	194	SetSigmaTau(3.0); // sigma_tau = 3 fm/c
	195	SetRmaxB(11.0); // fR = 11 fm
	196	SetYlMax(4.0); // fYmax = 4.0
	197	SetEtaRMax(1.1); // Umax = 1.1
	198	SetMomAsymmPar(0.0); // delta = 0.0
	199	SetCoordAsymmPar(0.0); // epsilon = 0.0
	200	SetFlagWeakDecay(0); // weak decay on (<0 off !!!)
	201	SetEtaType(1); // gaus distributed with fYmax dispersion (0 means boost invariant)
	202	SetGammaS(1.0); // gammaS = 1.0 (no strangeness canonical suppresion)
	203	SetPyquenNhsel(2); // hydro on, jets on, jet quenching on
	204	SetPyquenShad(1); // shadowing on (0 off)
	205	SetPyquenPtmin(7.0); // ptmin = 7.0 GeV/c
	206	SetPyquenT0(0.8); // T0 = 800 MeV
	207	SetPyquenTau0(0.1); // tau0 = 0.4 fm/c
	208	SetPyquenNf(0); // 0 flavours
	209	SetPyquenIenglu(0); // radiative and collisional energy loss
	210	SetPyquenIanglu(0); // small gluon angular distribution
	211	}
	212
7b7936e9	213	TObjArray* TUHKMgen::ImportParticles(const Option_t *)
b1c2e580	214	{
7b7936e9	215	// Function overloading the TGenerator::ImportParticles() member function.
	216	// The particles from the local particle list (fSecondariesList) are
	217	// forwarded to the TGenerator::fParticles
b1c2e580	218	fParticles->Clear();
b1c2e580	219	Int_t nump = 0;
b1c2e580	220	LPIT_t it,e;
b1c2e580	221
3fa37a65	222	for(it = fSecondariesList.begin(), e = fSecondariesList.end(); it != e; it++) {
b1c2e580	223	TVector3 pos(it->Pos().Vect());
	224	TVector3 mom(it->Mom().Vect());
	225	Float_t m1 = it->TableMass();
	226	Int_t im1 = it->GetMother();
	227	Int_t im2 = -1;
	228	Int_t id1 = -1;
	229	Int_t id2 = -1;
b1c2e580	230
b1c2e580	231	Int_t type = it->GetType(); // 0-hydro, 1-jets
	232
	233	if (im1> -1) {
3fa37a65	234	TParticle mother = (TParticle) (fParticles->UncheckedAt(im1+1));
b1c2e580	235	mother->SetLastDaughter(nump);
b1c2e580	236	if (mother->GetFirstDaughter()==-1)
3fa37a65	237	mother->SetFirstDaughter(nump+1);
b1c2e580	238	}
	239
	240	nump++;
	241
	242	TParticle* p = new TParticle(it->Encoding(), type, //pdg,stat
	243	im1, im2, id1, id2, //m1,m2,d1,d2
	244	mom[0], mom[1], mom[2], TMath::Sqrt(mom.Mag2() + m1 * m1), //px,py,pz,e
7b7936e9	245	pos[0]1.e-13, pos[1]1.e-13, pos[2]*1.e-13,
7b7936e9	246	it->T()*1.e-13/3e10); //x,y,z,t
b1c2e580	247
	248	p->SetUniqueID(nump);
	249	fParticles->Add(p);
	250	} //end for
	251
b1c2e580	252	fAllocator.FreeList(fSecondariesList);
3fa37a65	253
b1c2e580	254	return fParticles;
	255	}
	256
7b7936e9	257	Int_t TUHKMgen::ImportParticles(TClonesArray particles, const Option_t option)
b1c2e580	258	{
7b7936e9	259	// Function overloading the TGenerator::ImportParticles() member function.
	260	// The particles from the local particle list (fSecondariesList) are
	261	// forwarded to the TGenerator::fParticles
3fa37a65	262
	263
	264
b1c2e580	265	if(particles==0) return 0;
	266	TClonesArray &particlesR=*particles;
	267	particlesR.Clear();
b1c2e580	268
	269	Int_t numprim,numsec; numprim=numsec=0;
	270	Int_t nump = 0;
	271	LPIT_t it,e;
7b7936e9	272
3fa37a65	273	for(it = fSecondariesList.begin(), e = fSecondariesList.end(); it != e; it++) {
	274	TVector3 pos(it->Pos().Vect());
	275	TVector3 mom(it->Mom().Vect());
	276	Float_t m1 = it->TableMass();
	277	Int_t im1 = it->GetMother();
b1c2e580	278	Int_t im2 = -1;
	279	Int_t id1 = -1;
	280	Int_t id2 = -1;
b1c2e580	281
3fa37a65	282	Int_t type = it->GetType();
3fa37a65	283
b1c2e580	284	if (im1> -1) {
3fa37a65	285	// particle not a primary -> set the daughter indexes for the mother particle"<< endl;
3fa37a65	286	TParticle mother = (TParticle) (particlesR.UncheckedAt(im1));
b1c2e580	287	mother->SetLastDaughter(nump);
3fa37a65	288	if(mother->GetFirstDaughter()==-1)
3fa37a65	289	mother->SetFirstDaughter(nump);
b1c2e580	290	}
b1c2e580	291
3fa37a65	292
b1c2e580	293	new (particlesR[nump]) TParticle(it->Encoding(), type, //pdg,stat
	294	im1, im2, id1, id2, //m1,m2,d1,d2
	295	mom[0], mom[1], mom[2], TMath::Sqrt(mom.Mag2() + m1 * m1), //px,py,pz,e
7b7936e9	296	pos[0]1.e-13, pos[1]1.e-13, pos[2]*1.e-13,
	297	it->T()*1.e-13/3e10); //x,y,z,t
	298
b1c2e580	299	particlesR[nump]->SetUniqueID(nump);
3fa37a65	300	nump++;
b1c2e580	301	numsec++;
	302	}//end for
	303
3fa37a65	304	fSecondariesList.clear();
	305	printf("Scan and add prim %d sec %d and all %d particles\n",
	306	numprim,numsec,nump);
b1c2e580	307	return nump;
	308	}
	309
	310	//______________________________________________________________________________
	311	void TUHKMgen::Initialize()
	312	{
7b7936e9	313	// Function overloading the TGenerator::Initialize() member function.
	314	// The Monte-Carlo model is initialized (input parameters are transmitted,
	315	// particle list and decay channels are loaded, average multiplicities are calculated, etc.)
3fa37a65	316
	317	fInitialState = new InitialStateHydjet();
	318	SetAllParameters();
	319	fInitialState->LoadPDGInfo();
b1c2e580	320	// set the stable flags
3fa37a65	321	for(Int_t i=0; i<fStableFlagged; i++)
3fa37a65	322	fInitialState->SetPDGParticleStable(fStableFlagPDG[i], fStableFlagStatus[i]);
b1c2e580	323
3fa37a65	324	if(!fInitialState->MultIni())
7b7936e9	325	Error("TUHKMgen::Initialize()", "Bad status return from MultIni(). Check it out!! \n"); //
b1c2e580	326	}
b1c2e580	327
786056a2	328	void TUHKMgen::Print(const Option_t*) const
b1c2e580	329	{
	330	cout << "TUHKMgen::Print() method not implemented yet!!" << endl;
	331	}
	332
	333	void TUHKMgen::GenerateEvent()
	334	{
7b7936e9	335	// Member function overloading the TGenerator::GenerateEvent()
7b7936e9	336	// The HYDJET++ model is run and the particle lists (fSourceList and fSecondariesList) are filled
7b7936e9	337
3fa37a65	338	fInitialState->Initialize(fSecondariesList, fAllocator);
b1c2e580	339
3fa37a65	340	if(fSecondariesList.empty())Error("TUHKMgen::GenerateEvent()", "Source particle list empty after fireball initialization!! \n"); //
b1c2e580	341
b1c2e580	342	// Run the decays
7b7936e9	343
3fa37a65	344	if(fInitialState->RunDecays())
	345	fInitialState->Evolve(fSecondariesList, fAllocator, fInitialState->GetWeakDecayLimit());
	346
b1c2e580	347	}
	348
	349	void TUHKMgen::SetAllParameters() {
7b7936e9	350	// forward all model input parameters to the InitialStateHydjet object
	351	// which will handle all the Monte-Carlo simulation using HYDJET++ model
	352
b1c2e580	353	fInitialState->fParams.fSqrtS = fHydjetParams.fSqrtS;
	354	fInitialState->fParams.fAw = fHydjetParams.fAw;
	355	fInitialState->fParams.fIfb = fHydjetParams.fIfb ;
	356	fInitialState->fParams.fBfix = fHydjetParams.fBfix;
	357	fInitialState->fParams.fBmin = fHydjetParams.fBmin;
	358	fInitialState->fParams.fBmax = fHydjetParams.fBmax;
	359	fInitialState->fParams.fSeed = fHydjetParams.fSeed;
	360	fInitialState->fParams.fT = fHydjetParams.fT;
	361	fInitialState->fParams.fMuB = fHydjetParams.fMuB;
	362	fInitialState->fParams.fMuS = fHydjetParams.fMuS;
	363	fInitialState->fParams.fMuI3 = fHydjetParams.fMuI3;
	364	fInitialState->fParams.fThFO = fHydjetParams.fThFO;
	365	fInitialState->fParams.fMu_th_pip = fHydjetParams.fMu_th_pip;
	366
	367	fInitialState->fParams.fTau = fHydjetParams.fTau;
	368	fInitialState->fParams.fSigmaTau = fHydjetParams.fSigmaTau;
	369	fInitialState->fParams.fR = fHydjetParams.fR;
	370	fInitialState->fParams.fYlmax = fHydjetParams.fYlmax;
	371	fInitialState->fParams.fUmax = fHydjetParams.fUmax;
	372	fInitialState->fParams.fDelta = fHydjetParams.fDelta;
	373	fInitialState->fParams.fEpsilon = fHydjetParams.fEpsilon;
	374
b1c2e580	375	fInitialState->fParams.fWeakDecay = fHydjetParams.fWeakDecay;
3fa37a65	376	fInitialState->fParams.fDecay = 1; // run decays
b1c2e580	377	fInitialState->fParams.fCorrS = fHydjetParams.fCorrS;
	378	fInitialState->fParams.fEtaType = fHydjetParams.fEtaType;
	379
	380	fInitialState->fParams.fPtmin = fHydjetParams.fPtmin;
	381	fInitialState->fParams.fNhsel = fHydjetParams.fNhsel;
	382	fInitialState->fParams.fIshad = fHydjetParams.fIshad;
	383	fInitialState->fParams.fT0 = fHydjetParams.fT0;
	384	fInitialState->fParams.fTau0 = fHydjetParams.fTau0;
	385	fInitialState->fParams.fNf = fHydjetParams.fNf;
	386	fInitialState->fParams.fIenglu = fHydjetParams.fIenglu;
	387	fInitialState->fParams.fIanglu = fHydjetParams.fIanglu;
	388
	389	fInitialState->SetPDGParticleFilename(fParticleFilename);
	390	fInitialState->SetPDGDecayFilename(fDecayFilename);
7b7936e9	391	// fInitialState->SetUseCharmParticles(fUseCharmParticles);
	392	// fInitialState->SetWidthRange(fMinWidth, fMaxWidth);
	393	// fInitialState->SetMassRange(fMinMass, fMaxMass);
	394	// for(Int_t i=0; i<fStableFlagged; i++)
	395	// fInitialState->SetPDGParticleStable(fStableFlagPDG[i], fStableFlagStatus[i]);
	396	// cout << "TUHKMgen::SetAllParameters() OUT" << endl;
b1c2e580	397	}
b1c2e580	398