[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),
  fAllocator(),
  fSecondariesList(),
  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;  
*/

  strncpy(fParticleFilename, Form("%s/TUHKMgen/UHKM/particles.data", gSystem->Getenv("ALICE_ROOT")), 255);
  strncpy(fDecayFilename, Form("%s/TUHKMgen/UHKM/tabledecay.txt", gSystem->Getenv("ALICE_ROOT")), 255);
  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
  option = option;   // just to avoid the warning
  

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