]>
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; |
34 | ||
b1c2e580 | 35 | ClassImp(TUHKMgen) |
36 | ||
37 | TUHKMgen::TUHKMgen() : | |
38 | TGenerator("UHKM","UHKM"), | |
39 | fInitialState(0x0), | |
25c1936c | 40 | fAllocator(), |
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) |
48 | ||
25c1936c | 49 | // ParticleAllocator fAllocator; |
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); |
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) |
139 | delete fInitialState; | |
b1c2e580 | 140 | } |
141 | ||
142 | void TUHKMgen::SetAllParametersRHIC() | |
143 | { | |
7b7936e9 | 144 | // Set reasonable input parameters for 0-5% central Au+Au collisions |
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; |
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); |
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, |
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 | |
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(); |
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) |
291 | mother->SetFirstDaughter(nump); | |
7fecd07a | 292 | } else |
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++) |
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 | } |
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() |
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 | |
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 | } |