1 /////////////////////////////////////////////////////////////////////////////
2 // TUHKM is an interface to
4 // ( only the HYDJET++ part is implemented) //
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 //
16 ////////////////////////////////////////////////////////////////////////////
21 #include "TObjArray.h"
22 #include "TParticle.h"
26 #include "TClonesArray.h"
37 TUHKMgen::TUHKMgen() :
38 TGenerator("UHKM","UHKM"),
47 // default constructor setting reasonable defaults for initial parameters (central Pb+Pb at 5.5 TeV)
49 // ParticleAllocator fAllocator;
50 // List_t fSecondariesList;
52 // Set reasonable default values for LHC
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.;
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;
85 fHydjetParams.fIenglu=0;
86 fHydjetParams.fIanglu=0;
89 // Set reasonable default values for RHIC
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;
120 fHydjetParams.fIenglu=0;
121 fHydjetParams.fIanglu=0;
124 strncpy(fParticleFilename, Form("%s/TUHKMgen/UHKM/particles.data", gSystem->Getenv("ALICE_ROOT")), 255);
125 strncpy(fDecayFilename, Form("%s/TUHKMgen/UHKM/tabledecay.txt", gSystem->Getenv("ALICE_ROOT")), 255);
126 for(Int_t i=0; i<500; i++) {
127 fStableFlagPDG[i] = 0;
128 fStableFlagStatus[i] = kFALSE;
131 // cout << "TUHKMgen::TUHKMgen() OUT" << endl;
134 //______________________________________________________________________________
135 TUHKMgen::~TUHKMgen()
137 // destructor, deletes the InitialStateHydjet object
139 delete fInitialState;
142 void TUHKMgen::SetAllParametersRHIC()
144 // Set reasonable input parameters for 0-5% central Au+Au collisions
145 // at 200 GeV at RHIC
146 SetEcms(200.0); // RHIC top energy
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
179 void TUHKMgen::SetAllParametersLHC()
181 // Set reasonable input parameters for 0-5% Pb+Pb collisions at 5.5 TeV at LHC
182 SetEcms(5500.0); // LHC
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
215 TObjArray* TUHKMgen::ImportParticles(const Option_t *)
217 // Function overloading the TGenerator::ImportParticles() member function.
218 // The particles from the local particle list (fSecondariesList) are
219 // forwarded to the TGenerator::fParticles
224 for(it = fSecondariesList.begin(), e = fSecondariesList.end(); it != e; it++) {
225 TVector3 pos(it->Pos().Vect());
226 TVector3 mom(it->Mom().Vect());
227 Float_t m1 = it->TableMass();
228 Int_t im1 = it->GetMother();
233 Int_t type = it->GetType(); // 0-hydro, 1-jets
236 TParticle *mother = (TParticle*) (fParticles->UncheckedAt(im1+1));
237 mother->SetLastDaughter(nump);
238 if (mother->GetFirstDaughter()==-1)
239 mother->SetFirstDaughter(nump+1);
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
247 pos[0]*1.e-13, pos[1]*1.e-13, pos[2]*1.e-13,
248 it->T()*1.e-13/3e10); //x,y,z,t
250 p->SetUniqueID(nump);
254 fAllocator.FreeList(fSecondariesList);
259 Int_t TUHKMgen::ImportParticles(TClonesArray *particles, const Option_t* option)
261 // Function overloading the TGenerator::ImportParticles() member function.
262 // The particles from the local particle list (fSecondariesList) are
263 // forwarded to the TGenerator::fParticles
264 option = option; // just to avoid the warning
267 if(particles==0) return 0;
268 TClonesArray &particlesR=*particles;
271 Int_t numprim,numsec; numprim=numsec=0;
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();
284 Int_t type = it->GetType();
287 // particle not a primary -> set the daughter indexes for the mother particle"<< endl;
288 TParticle *mother = (TParticle*) (particlesR.UncheckedAt(im1));
289 mother->SetLastDaughter(nump);
290 if(mother->GetFirstDaughter()==-1)
291 mother->SetFirstDaughter(nump);
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
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
301 particlesR[nump]->SetUniqueID(nump);
306 fSecondariesList.clear();
307 printf("Scan and add prim %d sec %d and all %d particles\n",
308 numprim,numsec,nump);
312 //______________________________________________________________________________
313 void TUHKMgen::Initialize()
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.)
319 fInitialState = new InitialStateHydjet();
321 fInitialState->LoadPDGInfo();
322 // set the stable flags
323 for(Int_t i=0; i<fStableFlagged; i++)
324 fInitialState->SetPDGParticleStable(fStableFlagPDG[i], fStableFlagStatus[i]);
326 if(!fInitialState->MultIni())
327 Error("TUHKMgen::Initialize()", "Bad status return from MultIni(). Check it out!! \n"); //
330 void TUHKMgen::Print(const Option_t*) const
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;
340 void TUHKMgen::GenerateEvent()
342 // Member function overloading the TGenerator::GenerateEvent()
343 // The HYDJET++ model is run and the particle lists (fSourceList and fSecondariesList) are filled
345 fInitialState->Initialize(fSecondariesList, fAllocator);
347 if(fSecondariesList.empty())Error("TUHKMgen::GenerateEvent()", "Source particle list empty after fireball initialization!! \n"); //
351 if(fInitialState->RunDecays())
352 fInitialState->Evolve(fSecondariesList, fAllocator, fInitialState->GetWeakDecayLimit());
356 void TUHKMgen::SetAllParameters() {
357 // forward all model input parameters to the InitialStateHydjet object
358 // which will handle all the Monte-Carlo simulation using HYDJET++ model
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;
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;
382 fInitialState->fParams.fWeakDecay = fHydjetParams.fWeakDecay;
383 fInitialState->fParams.fDecay = 1; // run decays
384 fInitialState->fParams.fCorrS = fHydjetParams.fCorrS;
385 fInitialState->fParams.fEtaType = fHydjetParams.fEtaType;
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;
396 fInitialState->SetPDGParticleFilename(fParticleFilename);
397 fInitialState->SetPDGDecayFilename(fDecayFilename);
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;