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 // fUseCharmParticles(kFALSE),
53 // default constructor setting reasonable defaults for initial parameters (central Pb+Pb at 5.5 TeV)
55 // cout << "TUHKMgen::TUHKMgen() IN" << endl;
57 // Set reasonable default values for LHC
59 fHydjetParams.fSqrtS=5500; //LHC
60 fHydjetParams.fAw=207;//Au-Au
61 fHydjetParams.fIfb = 1; // flag of type of centrality generation (=0 is fixed by fBfix, not 0
62 //impact parameter is generated in each event between fBfmin
63 //and fBmax according with Glauber model (f-la 30)
64 fHydjetParams.fBfix=0.;
65 fHydjetParams.fBmin=0.;
66 fHydjetParams.fBmax=0.5; //0-5% centrality
67 fHydjetParams.fT = 0.170;
68 fHydjetParams.fMuB = 0.0;
69 fHydjetParams.fMuS = 0.0;
70 fHydjetParams.fMuI3 = 0.0;
71 fHydjetParams.fThFO = 0.130;
72 fHydjetParams.fMu_th_pip = 0.0;
73 fHydjetParams.fSeed=0;
74 fHydjetParams.fTau=10.;
75 fHydjetParams.fSigmaTau=3.;
77 fHydjetParams.fYlmax=4.0;
78 fHydjetParams.fUmax=1.1;
79 fHydjetParams.fDelta=0.;
80 fHydjetParams.fEpsilon=0.;
81 fHydjetParams.fWeakDecay=0; //>=0 on ,-1 off
82 fHydjetParams.fEtaType=1;//gaus
83 fHydjetParams.fCorrS=1.;
84 fHydjetParams.fNhsel=2;
85 fHydjetParams.fIshad=1;
86 fHydjetParams.fPtmin=7.0;
87 fHydjetParams.fT0=0.8;
88 fHydjetParams.fTau0=0.1;
90 fHydjetParams.fIenglu=0;
91 fHydjetParams.fIanglu=0;
94 // Set reasonable default values for RHIC
96 fHydjetParams.fSqrtS=200; //RHIC
97 fHydjetParams.fAw=197;//Au-Au
98 fHydjetParams.fIfb = 1;
99 fHydjetParams.fBfix=0.;
100 fHydjetParams.fBmin=0.;
101 fHydjetParams.fBmax=0.5; //0-5% centrality
102 fHydjetParams.fT = 0.165;
103 fHydjetParams.fMuB = 0.0285;
104 fHydjetParams.fMuS = 0.007;
105 fHydjetParams.fMuI3 = -0.001;
106 fHydjetParams.fThFO = 0.100;
107 fHydjetParams.fMu_th_pip = 0.053;
108 fHydjetParams.fSeed=0;
109 fHydjetParams.fTau=8.;
110 fHydjetParams.fSigmaTau=2.;
111 fHydjetParams.fR=10.;
112 fHydjetParams.fYlmax=3.3;
113 fHydjetParams.fUmax=1.1;
114 fHydjetParams.fDelta=0.;
115 fHydjetParams.fEpsilon=0.;
116 fHydjetParams.fWeakDecay=0;//>=0 on ,-1 off
117 fHydjetParams.fEtaType=1;//gaus
118 fHydjetParams.fCorrS=1.;
119 fHydjetParams.fNhsel=2;
120 fHydjetParams.fIshad=1;
121 fHydjetParams.fPtmin=3.4;
122 fHydjetParams.fT0=0.3;
123 fHydjetParams.fTau0=0.4;
125 fHydjetParams.fIenglu=0;
126 fHydjetParams.fIanglu=0;
129 strcpy(fParticleFilename, Form("%s/TUHKMgen/UHKM/particles.data", gSystem->Getenv("ALICE_ROOT")));
130 strcpy(fDecayFilename, Form("%s/TUHKMgen/UHKM/tabledecay.txt", gSystem->Getenv("ALICE_ROOT")));
131 for(Int_t i=0; i<500; i++) {
132 fStableFlagPDG[i] = 0;
133 fStableFlagStatus[i] = kFALSE;
136 // cout << "TUHKMgen::TUHKMgen() OUT" << endl;
139 //______________________________________________________________________________
140 TUHKMgen::~TUHKMgen()
142 // destructor, deletes the InitialStateHydjet object
143 // cout << "TUHKMgen::~TUHKMgen() IN" << endl;
145 delete fInitialState;
146 // cout << "TUHKMgen::~TUHKMgen() OUT" << endl;
149 void TUHKMgen::SetAllParametersRHIC()
151 // Set reasonable input parameters for 0-5% central Au+Au collisions
152 // at 200 GeV at RHIC
153 SetEcms(200.0); // RHIC top energy
159 SetChFrzTemperature(0.165); // T_ch = 165 MeV
160 SetMuB(0.0285); // mu_B = 28.5 MeV
161 SetMuS(0.007); // mu_S = 7 MeV
162 SetMuQ(-0.001); // mu_Q = -1 MeV
163 SetThFrzTemperature(0.100); // T_th = 100 MeV
164 SetMuPionThermal(0.053); // mu_th_pion = 53 MeV
165 SetSeed(0); // use UNIX time
166 SetTauB(8.0); // tau = 8 fm/c
167 SetSigmaTau(2.0); // sigma_tau = 2 fm/c
168 SetRmaxB(10.0); // fR = 10 fm
169 SetYlMax(3.3); // fYmax = 3.3
170 SetEtaRMax(1.1); // Umax = 1.1
171 SetMomAsymmPar(0.0); // delta = 0.0
172 SetCoordAsymmPar(0.0); // epsilon = 0.0
173 SetFlagWeakDecay(0.0); // weak decay on (<0 off !!!)
174 SetEtaType(1); // gaus distributed with fYmax dispersion (0 means boost invariant)
175 SetGammaS(1.0); // gammaS = 1.0 (no strangeness canonical suppresion)
176 SetPyquenNhsel(2); // hydro on, jets on, jet quenching on
177 SetPyquenShad(1); // shadowing on (0 off)
178 SetPyquenPtmin(3.4); // ptmin = 3.4 GeV/c
179 SetPyquenT0(0.3); // T0 = 300 MeV
180 SetPyquenTau0(0.4); // tau0 = 0.4 fm/c
181 SetPyquenNf(2); // 2 flavours
182 SetPyquenIenglu(0); // radiative and collisional energy loss
183 SetPyquenIanglu(0); // small gluon angular distribution
186 void TUHKMgen::SetAllParametersLHC()
188 // Set reasonable input parameters for 0-5% Pb+Pb collisions at 5.5 TeV at LHC
189 SetEcms(5500.0); // LHC
195 SetChFrzTemperature(0.170); // T_ch = 170 MeV
196 SetMuB(0.0); // mu_B = 0 MeV
197 SetMuS(0.0); // mu_S = 0 MeV
198 SetMuQ(0.0); // mu_Q = 0 MeV
199 SetThFrzTemperature(0.130); // T_th = 130 MeV
200 SetMuPionThermal(0.0); // mu_th_pion = 0 MeV
201 SetSeed(0); // use UNIX time
202 SetTauB(10.0); // tau = 10 fm/c
203 SetSigmaTau(3.0); // sigma_tau = 3 fm/c
204 SetRmaxB(11.0); // fR = 11 fm
205 SetYlMax(4.0); // fYmax = 4.0
206 SetEtaRMax(1.1); // Umax = 1.1
207 SetMomAsymmPar(0.0); // delta = 0.0
208 SetCoordAsymmPar(0.0); // epsilon = 0.0
209 SetFlagWeakDecay(0); // weak decay on (<0 off !!!)
210 SetEtaType(1); // gaus distributed with fYmax dispersion (0 means boost invariant)
211 SetGammaS(1.0); // gammaS = 1.0 (no strangeness canonical suppresion)
212 SetPyquenNhsel(2); // hydro on, jets on, jet quenching on
213 SetPyquenShad(1); // shadowing on (0 off)
214 SetPyquenPtmin(7.0); // ptmin = 7.0 GeV/c
215 SetPyquenT0(0.8); // T0 = 800 MeV
216 SetPyquenTau0(0.1); // tau0 = 0.4 fm/c
217 SetPyquenNf(0); // 0 flavours
218 SetPyquenIenglu(0); // radiative and collisional energy loss
219 SetPyquenIanglu(0); // small gluon angular distribution
222 TObjArray* TUHKMgen::ImportParticles(const Option_t *)
224 // Function overloading the TGenerator::ImportParticles() member function.
225 // The particles from the local particle list (fSecondariesList) are
226 // forwarded to the TGenerator::fParticles
227 // cout << "TObjArray* TUHKMgen::ImportParticles(Option_t) IN" << endl;
229 // cout << "TObjArray* TUHKMgen::ImportParticles(Option_t): UHKM stack contains " << fNPsec << " particles." << endl;
233 for(it = fSecondariesList.begin(), e = fSecondariesList.end(); it != e; ++it) {
234 TVector3 pos(it->Pos().Vect());
235 TVector3 mom(it->Mom().Vect());
236 Float_t m1 = it->TableMass();
237 Int_t im1 = it->GetMother();
242 // Int_t nd = it->GetNDaughters();
244 Int_t type = it->GetType(); // 0-hydro, 1-jets
247 TParticle *mother = (TParticle*) (fParticles->UncheckedAt(im1));
248 mother->SetLastDaughter(nump);
249 if (mother->GetFirstDaughter()==-1)
250 mother->SetFirstDaughter(nump);
255 TParticle* p = new TParticle(it->Encoding(), type, //pdg,stat
256 im1, im2, id1, id2, //m1,m2,d1,d2
257 mom[0], mom[1], mom[2], TMath::Sqrt(mom.Mag2() + m1 * m1), //px,py,pz,e
258 pos[0]*1.e-13, pos[1]*1.e-13, pos[2]*1.e-13,
259 it->T()*1.e-13/3e10); //x,y,z,t
261 p->SetUniqueID(nump);
265 fAllocator.FreeList(fSourceList);
266 fAllocator.FreeList(fSecondariesList);
267 // cout << "TObjArray* TUHKMgen::ImportParticles(Option_t) OUT" << endl;
271 Int_t TUHKMgen::ImportParticles(TClonesArray *particles, const Option_t* option)
273 // Function overloading the TGenerator::ImportParticles() member function.
274 // The particles from the local particle list (fSecondariesList) are
275 // forwarded to the TGenerator::fParticles
276 // cout << "Int_t TUHKMgen::ImportParticles(TClonesArray*, Option_t) IN" << endl;
277 if(particles==0) return 0;
278 TClonesArray &particlesR=*particles;
282 Int_t numprim,numsec; numprim=numsec=0;
286 cout << "TUHKMgen::ImportParticles() option(All or Sec) = " << option << endl;
287 for(it = fSecondariesList.begin(), e = fSecondariesList.end(); it != e; ++it) {
288 //!!! for(Int_t pp=0;pp<100;pp++) {
289 // cout << "TUHKMgen::ImportParticles() import particle pdg(" << it->Encoding() << ")" << endl;
290 TVector3 pos(it->Pos().Vect()); //
291 //!!! TVector3 pos(0.0, 0.1, 0.2);
292 TVector3 mom(it->Mom().Vect()); //
293 //!!! TVector3 mom(-0.1, 0.1, 0.2);
294 Float_t m1 = it->TableMass(); //
295 //!!! Float_t m1 = 0.139;
296 Int_t im1 = it->GetMother(); //
297 //!!! Int_t im1 = (pp<50? -1 : pp-50);
303 // nd = it->GetNDaughters();
305 Int_t type = it->GetType(); // // 0-hydro, 1-jets
306 //!!! Int_t type = 0;
309 TParticle *mother = (TParticle*) (particlesR.UncheckedAt(im1+1));
310 mother->SetLastDaughter(nump);
311 if (mother->GetFirstDaughter()==-1)
312 mother->SetFirstDaughter(nump);
317 new (particlesR[nump]) TParticle(it->Encoding(), type, //pdg,stat
318 im1, im2, id1, id2, //m1,m2,d1,d2
319 mom[0], mom[1], mom[2], TMath::Sqrt(mom.Mag2() + m1 * m1), //px,py,pz,e
320 pos[0]*1.e-13, pos[1]*1.e-13, pos[2]*1.e-13,
321 it->T()*1.e-13/3e10); //x,y,z,t
323 /*!!! new (particlesR[nump]) TParticle(211, type, //pdg,stat
324 im1, im2, id1, id2, //m1,m2,d1,d2
325 mom[0], mom[1], mom[2], TMath::Sqrt(mom.Mag2() + m1 * m1), //px,py,pz,e
326 pos[0]*1.e-13, pos[1]*1.e-13, pos[2]*1.e-13,
327 0.0*1.e-13/3e10); //x,y,z,t
330 particlesR[nump]->SetUniqueID(nump);
334 fAllocator.FreeList(fSourceList); //
335 fAllocator.FreeList(fSecondariesList); //
336 // printf("Scan and add prim %d sec %d and all %d particles\n",
337 // numprim,numsec,nump);
338 // cout << "Int_t TUHKMgen::ImportParticles(TClonesArray*, Option_t) OUT" << endl;
342 //______________________________________________________________________________
343 void TUHKMgen::Initialize()
345 // Function overloading the TGenerator::Initialize() member function.
346 // The Monte-Carlo model is initialized (input parameters are transmitted,
347 // particle list and decay channels are loaded, average multiplicities are calculated, etc.)
348 // cout << "TUHKMgen::Initialize() IN" << endl;
349 fInitialState = new InitialStateHydjet(); //
350 SetAllParameters(); //
351 fInitialState->LoadPDGInfo(); //
352 // set the stable flags
353 for(Int_t i=0; i<fStableFlagged; i++) //
354 fInitialState->SetPDGParticleStable(fStableFlagPDG[i], fStableFlagStatus[i]); //
356 if(!fInitialState->MultIni()) //
357 Error("TUHKMgen::Initialize()", "Bad status return from MultIni(). Check it out!! \n"); //
358 // cout << "TUHKMgen::Initialize() OUT" << endl;
361 void TUHKMgen::Print(const Option_t*) const
363 cout << "TUHKMgen::Print() method not implemented yet!!" << endl;
366 void TUHKMgen::GenerateEvent()
368 // Member function overloading the TGenerator::GenerateEvent()
369 // The HYDJET++ model is run and the particle lists (fSourceList and fSecondariesList) are filled
370 // cout << "TUHKMgen::GenerateEvent() IN" << endl;
371 // cout << "TUHKMgen::GenerateEvent() IN fSourceList " <<endl;
372 // cout<< " fInitialStrate "<<fInitialState<<endl; //
373 //" fSourceList "<< fSourceList<<" fAllocator "<< fAllocator <<endl;
375 // Generate the initial particles
376 fInitialState->Initialize(fSourceList, fAllocator); //
377 // cout << "TUHKMgen::fInitialState->Initialize" << endl;
379 if(fSourceList.empty()) //
380 Error("TUHKMgen::GenerateEvent()", "Source particle list empty after fireball initialization!! \n"); //
384 // cout << "after Initilize: get_time, weak_decay_limit" <<fInitialState->GetTime()<<fInitialState->GetWeakDecayLimit()<< endl; //
385 if(fInitialState->GetTime() >= 0.) //
386 fInitialState->Evolve(fSourceList, fSecondariesList, fAllocator, fInitialState->GetWeakDecayLimit()); //
387 // cout << "TUHKMgen::GenerateEvent() OUT" << endl;
390 void TUHKMgen::SetAllParameters() {
391 // forward all model input parameters to the InitialStateHydjet object
392 // which will handle all the Monte-Carlo simulation using HYDJET++ model
394 // cout << "TUHKMgen::SetAllParameters() IN" << endl;
396 fInitialState->fParams.fSqrtS = fHydjetParams.fSqrtS;
397 fInitialState->fParams.fAw = fHydjetParams.fAw;
398 fInitialState->fParams.fIfb = fHydjetParams.fIfb ;
399 fInitialState->fParams.fBfix = fHydjetParams.fBfix;
400 fInitialState->fParams.fBmin = fHydjetParams.fBmin;
401 fInitialState->fParams.fBmax = fHydjetParams.fBmax;
402 fInitialState->fParams.fSeed = fHydjetParams.fSeed;
403 fInitialState->fParams.fT = fHydjetParams.fT;
404 fInitialState->fParams.fMuB = fHydjetParams.fMuB;
405 fInitialState->fParams.fMuS = fHydjetParams.fMuS;
406 fInitialState->fParams.fMuI3 = fHydjetParams.fMuI3;
407 fInitialState->fParams.fThFO = fHydjetParams.fThFO;
408 fInitialState->fParams.fMu_th_pip = fHydjetParams.fMu_th_pip;
410 fInitialState->fParams.fTau = fHydjetParams.fTau;
411 fInitialState->fParams.fSigmaTau = fHydjetParams.fSigmaTau;
412 fInitialState->fParams.fR = fHydjetParams.fR;
413 fInitialState->fParams.fYlmax = fHydjetParams.fYlmax;
414 fInitialState->fParams.fUmax = fHydjetParams.fUmax;
415 fInitialState->fParams.fDelta = fHydjetParams.fDelta;
416 fInitialState->fParams.fEpsilon = fHydjetParams.fEpsilon;
418 // fInitialState->fParams.fTime = fHydjetParams.fTime;
419 fInitialState->fParams.fWeakDecay = fHydjetParams.fWeakDecay;
420 fInitialState->fParams.fCorrS = fHydjetParams.fCorrS;
421 fInitialState->fParams.fEtaType = fHydjetParams.fEtaType;
423 fInitialState->fParams.fPtmin = fHydjetParams.fPtmin;
424 fInitialState->fParams.fNhsel = fHydjetParams.fNhsel;
425 fInitialState->fParams.fIshad = fHydjetParams.fIshad;
426 fInitialState->fParams.fT0 = fHydjetParams.fT0;
427 fInitialState->fParams.fTau0 = fHydjetParams.fTau0;
428 fInitialState->fParams.fNf = fHydjetParams.fNf;
429 fInitialState->fParams.fIenglu = fHydjetParams.fIenglu;
430 fInitialState->fParams.fIanglu = fHydjetParams.fIanglu;
432 fInitialState->SetPDGParticleFilename(fParticleFilename);
433 fInitialState->SetPDGDecayFilename(fDecayFilename);
434 // fInitialState->SetUseCharmParticles(fUseCharmParticles);
435 // fInitialState->SetWidthRange(fMinWidth, fMaxWidth);
436 // fInitialState->SetMassRange(fMinMass, fMaxMass);
437 // for(Int_t i=0; i<fStableFlagged; i++)
438 // fInitialState->SetPDGParticleStable(fStableFlagPDG[i], fStableFlagStatus[i]);
439 // cout << "TUHKMgen::SetAllParameters() OUT" << endl;