Updated versions (Adrian Sevcenco)
[u/mrichter/AliRoot.git] / TUHKMgen / AliGenUHKM.cxx
CommitLineData
b1c2e580 1/////////////////////////////////////////////////////////////////////////////
2// Generator using UHKM 3.0 as an external generator. //
3// ( only the HYDJET++ part is implemented for a moment) //
4// temporary link: //
5// http://lav01.sinp.msu.ru/~igor/hydjet++/hydjet++.txt //
6// The main UHKM options are accessable through this interface. //
7// Uses the TUHKMgen implementation of TGenerator. //
8// Author of the first implementation: Sergey Zaporozhets //
9// (zaporozh@sunhe.jinr.ru) //
10// Futhers modifications were made by //
11// Ionut Cristian Arsene (i.c.arsene@fys.uio.no) //
12// & Malinina Liudmila(malinina@lav01.sinp.msu.ru) using as an example //
13// AliGenTherminator.cxx created by Adam Kisiel //
14// //
15////////////////////////////////////////////////////////////////////////////
16
7b7936e9 17#include "TUHKMgen.h"
b1c2e580 18#ifndef DATABASE_PDG
19#include "DatabasePDG.h"
20#endif
21#ifndef PARTICLE_PDG
22#include "ParticlePDG.h"
23#endif
24#include <TLorentzVector.h>
25#include <TPDGCode.h>
26#include <TParticle.h>
27#include <TClonesArray.h>
28#include <TMCProcess.h>
29#include "TDatabasePDG.h"
30#include "TSystem.h"
31
32#include "AliGenUHKM.h"
33#include "AliRun.h"
34#include "AliConst.h"
35#include "AliDecayer.h"
36#include "AliGenEventHeader.h"
37#include "AliGenHijingEventHeader.h"
38#include "AliLog.h"
39
40#include <iostream>
41#include <string>
42using namespace std;
43
7b7936e9 44//class TClonesArray;
45//enum TMCProcess;
46//class AliRun;
47//class TSystem;
48//class AliDecayer;
49
b1c2e580 50ClassImp(AliGenUHKM)
51
52//_______________________________________
53AliGenUHKM::AliGenUHKM()
54 :AliGenMC(),
55 fTrials(0),
56 fUHKMgen(0),
786056a2 57 fHydjetParams(),
7b7936e9 58 fStableFlagged(0)
b1c2e580 59{
7b7936e9 60 // Default constructor setting up default reasonable parameter values
61 // for central Pb+Pb collisions at 5.5TeV
62 // cout << "AliGenUHKM::AliGenUHKM() IN" << endl;
b1c2e580 63
64 // LHC
65 fHydjetParams.fSqrtS=5500; //LHC
66 fHydjetParams.fAw=207;//Pb-Pb
67 fHydjetParams.fBmin=0.;
68 fHydjetParams.fBmax=0.5; //0-5% centrality
69 fHydjetParams.fT = 0.170;
70 fHydjetParams.fMuB = 0.0;
71 fHydjetParams.fMuS = 0.0;
72 fHydjetParams.fMuI3 = 0.0;
73 fHydjetParams.fThFO = 0.130;
74 fHydjetParams.fMu_th_pip = 0.0;
75 fHydjetParams.fSeed=0;
76 fHydjetParams.fTau=10.;
77 fHydjetParams.fSigmaTau=3.;
78 fHydjetParams.fR=11.;
79 fHydjetParams.fYlmax=4.0;
80 fHydjetParams.fUmax=1.1;
81 fHydjetParams.fDelta=0.;
82 fHydjetParams.fEpsilon=0.;
83 fHydjetParams.fWeakDecay=0; //>=0 on ,-1 off
84 fHydjetParams.fEtaType=1;//gaus
85 fHydjetParams.fCorrS=1.;
86 fHydjetParams.fNhsel=2;
87 fHydjetParams.fIshad=1;
88 fHydjetParams.fPtmin=7.0;
89 fHydjetParams.fT0=0.8;
90 fHydjetParams.fTau0=0.1;
91 fHydjetParams.fNf=0;
92 fHydjetParams.fIenglu=0;
93 fHydjetParams.fIanglu=0;
94
95
96/* RHIC
97 fHydjetParams.fSqrtS=200; //RHIC
98 fHydjetParams.fAw=197;//Au-Au
99 fHydjetParams.fBmin=0.;
100 fHydjetParams.fBmax=0.5; //0-5% centrality
101 fHydjetParams.fT = 0.165;
102 fHydjetParams.fMuB = 0.0285;
103 fHydjetParams.fMuS = 0.007;
104 fHydjetParams.fMuI3 = -0.001;
105 fHydjetParams.fThFO = 0.100;
106 fHydjetParams.fMu_th_pip = 0.053;
107 fHydjetParams.fSeed=0;
108 fHydjetParams.fTau=8.;
109 fHydjetParams.fSigmaTau=2.;
110 fHydjetParams.fR=10.;
111 fHydjetParams.fYlmax=3.3;
112 fHydjetParams.fUmax=1.1;
113 fHydjetParams.fDelta=0.;
114 fHydjetParams.fEpsilon=0.;
115 fHydjetParams.fWeakDecay=0; //>=0 on ,-1 off
116 fHydjetParams.fEtaType=1;//gaus
117 fHydjetParams.fCorrS=1.;
118 fHydjetParams.fNhsel=2;
119 fHydjetParams.fIshad=0;
120 fHydjetParams.fPtmin=3.4;
121 fHydjetParams.fT0=0.3;
122 fHydjetParams.fTau0=0.4;
123 fHydjetParams.fNf=2;
124 fHydjetParams.fIenglu=0;
125 fHydjetParams.fIanglu=0;
126*/
127 strcpy(fParticleFilename, Form("%s/TUHKMgen/UHKM/particles.data", gSystem->Getenv("ALICE_ROOT")));
128 strcpy(fDecayFilename, Form("%s/TUHKMgen/UHKM/tabledecay.txt", gSystem->Getenv("ALICE_ROOT")));
129 for(Int_t i=0; i<500; i++) {
130 fStableFlagPDG[i] = 0;
131 fStableFlagStatus[i] = kFALSE;
132 }
133 fStableFlagged = 0;
134
7b7936e9 135 // cout << "AliGenUHKM::AliGenUHKM() OUT" << endl;
b1c2e580 136}
137
138//_______________________________________
139AliGenUHKM::AliGenUHKM(Int_t npart)
140 :AliGenMC(npart),
141 fTrials(0),
142 fUHKMgen(0),
786056a2 143 fHydjetParams(),
7b7936e9 144 fStableFlagged(0)
b1c2e580 145{
7b7936e9 146 // Constructor specifying the size of the particle table
147 // and setting up default reasonable parameter values
148 // for central Pb+Pb collisions at 5.5TeV
149 // cout << "AliGenUHKM::AliGenUHKM(Int_t) IN" << endl;
b1c2e580 150 fName = "UHKM";
151 fTitle= "Particle Generator using UHKM 3.0";
152
b1c2e580 153 fNprimaries = 0;
154
155 //LHC
156 fHydjetParams.fSqrtS=5500; //LHC
157 fHydjetParams.fAw=207;//Pb-Pb
158 fHydjetParams.fBmin=0.;
159 fHydjetParams.fBmax=0.5; //0-5% centrality
160 fHydjetParams.fT = 0.170;
161 fHydjetParams.fMuB = 0.0;
162 fHydjetParams.fMuS = 0.0;
163 fHydjetParams.fMuI3 = 0.0;
164 fHydjetParams.fThFO = 0.130;
165 fHydjetParams.fMu_th_pip = 0.0;
166 fHydjetParams.fSeed=0;
167 fHydjetParams.fTau=10.;
168 fHydjetParams.fSigmaTau=3.;
169 fHydjetParams.fR=11.;
170 fHydjetParams.fYlmax=4.0;
171 fHydjetParams.fUmax=1.1;
172 fHydjetParams.fDelta=0.;
173 fHydjetParams.fEpsilon=0.;
174 fHydjetParams.fWeakDecay=0; //>=0 on ,-1 off
175 fHydjetParams.fEtaType=1;//gaus
176 fHydjetParams.fCorrS=1.;
177 fHydjetParams.fNhsel=2;
178 fHydjetParams.fIshad=1;
179 fHydjetParams.fPtmin=7.0;
180 fHydjetParams.fT0=0.8;
181 fHydjetParams.fTau0=0.1;
182 fHydjetParams.fNf=0;
183 fHydjetParams.fIenglu=0;
184 fHydjetParams.fIanglu=0;
185
186/*RHIC
187 fHydjetParams.fSqrtS=200; //RHIC
188 fHydjetParams.fAw=197;//Au-Au
189 fHydjetParams.fBmin=0.;
190 fHydjetParams.fBmax=0.5; //0-5% centrality
191 fHydjetParams.fT = 0.165;
192 fHydjetParams.fMuB = 0.0285;
193 fHydjetParams.fMuS = 0.007;
194 fHydjetParams.fMuI3 = -0.001;
195 fHydjetParams.fThFO = 0.100;
196 fHydjetParams.fMu_th_pip = 0.053;
197 fHydjetParams.fSeed=0;
198 fHydjetParams.fTau=8.;
199 fHydjetParams.fSigmaTau=2.;
200 fHydjetParams.fR=10.;
201 fHydjetParams.fYlmax=3.3;
202 fHydjetParams.fUmax=1.1;
203 fHydjetParams.fDelta=0.;
204 fHydjetParams.fEpsilon=0.;
205 fHydjetParams.fWeakDecay=0;//>=0 on ,-1 off
206 fHydjetParams.fEtaType=1;//gaus
207 fHydjetParams.fCorrS=1.;
208 fHydjetParams.fNhsel=2;
209 fHydjetParams.fIshad=1;
210 fHydjetParams.fPtmin=3.4;
211 fHydjetParams.fT0=0.3;
212 fHydjetParams.fTau0=0.4;
213 fHydjetParams.fNf=2;
214 fHydjetParams.fIenglu=0;
215 fHydjetParams.fIanglu=0;
216*/
217
218 strcpy(fParticleFilename, Form("%s/TUHKMgen/UHKM/particles.data", gSystem->Getenv("ALICE_ROOT")));
219 strcpy(fDecayFilename, Form("%s/TUHKMgen/UHKM/tabledecay.txt", gSystem->Getenv("ALICE_ROOT")));
220 for(Int_t i=0; i<500; i++) {
221 fStableFlagPDG[i] = 0;
222 fStableFlagStatus[i] = kFALSE;
223 }
224 fStableFlagged = 0;
225
7b7936e9 226 // cout << "AliGenUHKM::AliGenUHKM(Int_t) OUT" << endl;
b1c2e580 227}
228
229//__________________________________________
230AliGenUHKM::~AliGenUHKM()
7b7936e9 231{
232 // Destructor, do nothing
233}
b1c2e580 234
235void AliGenUHKM::SetAllParametersRHIC()
236{
7b7936e9 237 // Set reasonable default parameters for 0-5% central Au+Au collisions
238 // at 200 GeV at RHIC
b1c2e580 239 SetEcms(200.0); // RHIC top energy
240 SetAw(197); // Au+Au
241 SetBmin(0.0); // 0%
242 SetBmax(0.5); // 5%
243 SetChFrzTemperature(0.165); // T_ch = 165 MeV
244 SetMuB(0.0285); // mu_B = 28.5 MeV
245 SetMuS(0.007); // mu_S = 7 MeV
246 SetMuQ(-0.001); // mu_Q = -1 MeV
247 SetThFrzTemperature(0.100); // T_th = 100 MeV
248 SetMuPionThermal(0.053); // mu_th_pion = 53 MeV
249 SetSeed(0); // use UNIX time
250 SetTauB(8.0); // tau = 8 fm/c
251 SetSigmaTau(2.0); // sigma_tau = 2 fm/c
252 SetRmaxB(10.0); // fR = 10 fm
253 SetYlMax(3.3); // fYmax = 3.3
254 SetEtaRMax(1.1); // Umax = 1.1
255 SetMomAsymmPar(0.0); // delta = 0.0
256 SetCoordAsymmPar(0.0); // epsilon = 0.0
7b7936e9 257 // SetFlagWeakDecay(0); // weak decay on (<0 off !!!)
b1c2e580 258 SetEtaType(1); // gaus distributed with fYmax dispersion (0 means boost invariant)
259 SetGammaS(1.0); // gammaS = 1.0 (no strangeness canonical suppresion)
260 SetPyquenNhsel(2); // hydro on, jets on, jet quenching on
261 SetPyquenShad(1); // shadowing on (0 off)
262 SetPyquenPtmin(3.4); // ptmin = 3.4 GeV/c
263 SetPyquenT0(0.3); // T0 = 300 MeV
264 SetPyquenTau0(0.4); // tau0 = 0.4 fm/c
265 SetPyquenNf(2); // 2 flavours
266 SetPyquenIenglu(0); // radiative and collisional energy loss
267 SetPyquenIanglu(0); // small gluon angular distribution
268}
269
270void AliGenUHKM::SetAllParametersLHC()
271{
7b7936e9 272 // Set reasonable default parameters for 0-5% central Pb+Pb collisions
273 // at 5.5 TeV at LHC
b1c2e580 274 SetEcms(5500.0); // LHC
275 SetAw(207); // Pb+Pb
276 SetBmin(0.0); // 0%
277 SetBmax(0.5); // 5%
278 SetChFrzTemperature(0.170); // T_ch = 170 MeV
279 SetMuB(0.0); // mu_B = 0 MeV
280 SetMuS(0.0); // mu_S = 0 MeV
281 SetMuQ(0.0); // mu_Q = 0 MeV
282 SetThFrzTemperature(0.130); // T_th = 130 MeV
283 SetMuPionThermal(0.0); // mu_th_pion = 0 MeV
284 SetSeed(0); // use UNIX time
285 SetTauB(10.0); // tau = 10 fm/c
286 SetSigmaTau(3.0); // sigma_tau = 3 fm/c
287 SetRmaxB(11.0); // fR = 11 fm
288 SetYlMax(4.0); // fYmax = 4.0
289 SetEtaRMax(1.1); // Umax = 1.1
290 SetMomAsymmPar(0.0); // delta = 0.0
291 SetCoordAsymmPar(0.0); // epsilon = 0.0
7b7936e9 292 // SetFlagWeakDecay(0); // weak decay on (<0 off !!!)
b1c2e580 293 SetEtaType(1); // gaus distributed with fYmax dispersion (0 means boost invariant)
294 SetGammaS(1.0); // gammaS = 1.0 (no strangeness canonical suppresion)
295 SetPyquenNhsel(2); // hydro on, jets on, jet quenching on
296 SetPyquenShad(1); // shadowing on (0 off)
297 SetPyquenPtmin(7.0); // ptmin = 7.0 GeV/c
298 SetPyquenT0(0.8); // T0 = 800 MeV
299 SetPyquenTau0(0.1); // tau0 = 0.4 fm/c
300 SetPyquenNf(0); // 0 flavours
301 SetPyquenIenglu(0); // radiative and collisional energy loss
302 SetPyquenIanglu(0); // small gluon angular distribution
303}
304
305//_________________________________________
306void AliGenUHKM::Init()
307{
7b7936e9 308 // Initialization of the TGenerator::TUHKMgen interface object.
309 // Model input parameters are transmited to the TUHKMgen object which forwards them
310 // further to the model.
311 // HYDJET++ is initialized (average multiplicities are calculated, particle species definitions and decay
312 // channels are loaded, etc.)
313 // cout << "AliGenUHKM::Init() IN" << endl;
b1c2e580 314
315 SetMC(new TUHKMgen());
316 fUHKMgen = (TUHKMgen*) fMCEvGen;
317 SetAllParameters();
7b7936e9 318
b1c2e580 319 AliGenMC::Init();
320
321 fUHKMgen->Initialize();
7b7936e9 322 CheckPDGTable(); //
b1c2e580 323
7b7936e9 324 fUHKMgen->Print(); //
325 // cout << "AliGenUHKM::Init() OUT" << endl;
b1c2e580 326}
327
b1c2e580 328//________________________________________
329void AliGenUHKM::Generate()
330{
7b7936e9 331 // Generate one HYDJET++ event, get the output and push particles further
332 // to AliRoot's stack
333 // cout << "AliGenUHKM::Generate() IN" << endl;
b1c2e580 334 Float_t polar[3] = {0,0,0};
335 Float_t origin[3] = {0,0,0};
336 Float_t origin0[3] = {0,0,0};
337 Float_t v[3];
338 Float_t mass, energy;
339
340 Vertex();
341 for(Int_t j=0; j<3; j++) origin0[j] = fVertex[j];
342
343 fTrials = 0;
344 // cout << "AliGenUHKM::Generate() fTrials = " << fTrials << endl;
345
346 Int_t nt = 0;
347
348 fUHKMgen->GenerateEvent();
349 fTrials++;
350
351 fUHKMgen->ImportParticles(&fParticles,"All");
352
353 Int_t np = fParticles.GetEntriesFast();
7b7936e9 354 // cout << "AliGenUHKM::Generate() GetEntries " <<np<< endl;
b1c2e580 355
356
357 Int_t* idsOnStack = new Int_t[np];
358 Int_t* newPos = new Int_t[np];
359 for(Int_t i=0; i<np; i++) newPos[i] = i;
360
361 //_________ Loop for particle selection
362 // for(Int_t i=1; i<np; i++) {
363 for(Int_t i=1; i<np; i++) {
364 // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
365 // the particle indexes are 0 based but fParticles is a 1 based array
366 // -1 is the trivial code (when it does not exist)
367 // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
368 TParticle *iparticle = (TParticle*)fParticles.At(i);
7b7936e9 369 // cout << "AliGenUHKM::Generate() particle #" << i << " in fParticles *********************"<< endl;
b1c2e580 370
371 Int_t kf = iparticle->GetPdgCode();
7b7936e9 372 // if(kf==130 || kf==-130)
373 // cout << "AliGenUHKM::Generate() PDG = " << kf << endl;
b1c2e580 374
375 Bool_t hasMother = (iparticle->GetFirstMother() >= 0);
376
7b7936e9 377 // cout << "AliGenUHKM::Generate() mother index in fParticles = "
378 // << (iparticle->GetFirstMother()==-1 ? -1 : iparticle->GetFirstMother()+1)
379 // << " ; hasMother = " << hasMother << endl;
380 // cout << "AliGenUHKM::Generate() type (0-hydro;1-jet) = " << iparticle->GetStatusCode() << endl;
b1c2e580 381 Bool_t hasDaughter = (iparticle->GetNDaughters() > 0);
382
7b7936e9 383 // cout << "AliGenUHKM::Generate() n.daughters = " << iparticle->GetNDaughters()
b1c2e580 384 //<< " ; hasDaughter = " << hasDaughter << endl;
385
386
387 if(hasDaughter) {
388 // cout << "AliGenUHKM::Generate() decayed particle (not trackable)" << endl;
389 // This particle has decayed
390 // It will not be tracked
7b7936e9 391 // Add it only once with coordinates not
b1c2e580 392 // smeared with primary vertex position
393 Float_t p[3] = {p[0] = iparticle->Px(),
394 p[1] = iparticle->Py(),
395 p[2] = iparticle->Pz()};
396 mass = TDatabasePDG::Instance()->GetParticle(kf)->Mass();
397 energy = sqrt(mass*mass + p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
398 v[0] = iparticle->Vx();
399 v[1] = iparticle->Vy();
400 v[2] = iparticle->Vz();
401 Float_t time = iparticle->T();
402
7b7936e9 403 // if(kf==130 || kf==-130) cout << "type = " << iparticle->GetStatusCode() << endl; //j1/h0
b1c2e580 404
405 Int_t imo = -1;
406 if(hasMother) {
407 imo = iparticle->GetFirstMother(); //index of mother particle in fParticles
408 } // if has mother
7b7936e9 409 Bool_t trackFlag = kFALSE; // tFlag is kFALSE --> do not track the particle
b1c2e580 410
7b7936e9 411 // printf("Pushing Track %d with status %d mother %d\n", kf, trackFlag, imo>=0?idsOnStack[imo]:imo);
b1c2e580 412 PushTrack(trackFlag, (imo>=0 ? idsOnStack[imo+1] : imo), kf,
413 p[0], p[1], p[2], energy,
414 v[0], v[1], v[2], time,
415 polar[0], polar[1], polar[2],
416 (hasMother ? kPDecay : kPNoProcess), nt);
7b7936e9 417 // cout << "AliGenUHKM::Generate() pushed on stack with stack index = " << nt
418 // << "; mother index on stack = " << (imo>=0 ? idsOnStack[imo+1] : imo) << endl;
b1c2e580 419 idsOnStack[i] = nt;
420 fNprimaries++;
421 KeepTrack(nt);
422 }
423 else {
7b7936e9 424 // cout << "AliGenUHKM::Generate() final particle --> push it twice on the stack" << endl;
b1c2e580 425 // This is a final state particle
426 // It will be tracked
427 // Add it TWICE to the stack !!!
428 // First time with event-wide coordinates (for femtoscopy) -
429 // this one will not be tracked
7b7936e9 430 // Second time with event-wide c0ordinates and vertex smearing
b1c2e580 431 // this one will be tracked
432 Float_t p[3] = {p[0] = iparticle->Px(),
433 p[1] = iparticle->Py(),
434 p[2] = iparticle->Pz()};
435 mass = TDatabasePDG::Instance()->GetParticle(kf)->Mass();
436 energy = sqrt(mass*mass + p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
437 v[0] = iparticle->Vx();
438 v[1] = iparticle->Vy();
439 v[2] = iparticle->Vz();
440
7b7936e9 441 Int_t type = iparticle->GetStatusCode(); // 1-from jet / 0-from hydro
442 //if(kf==130 || kf==-130) cout << "type = " << iparticle->GetStatusCode() << endl; //j1/h0
b1c2e580 443 Int_t coeffT=1;
7b7936e9 444 if(type==1) coeffT=-1; //to separate particles from jets
b1c2e580 445
446 Int_t imo = -1;
447
448 if(hasMother) {
449 imo = iparticle->GetFirstMother();
450 } // if has mother
7b7936e9 451 Bool_t trackFlag = kFALSE; // tFlag = kFALSE --> do not track this one, its for femtoscopy
b1c2e580 452 PushTrack(trackFlag, (imo>=0 ? idsOnStack[imo+1] : imo), kf,
453 p[0], p[1], p[2], energy,
454 v[0], v[1], v[2], (iparticle->T())*coeffT,
455 polar[0], polar[1], polar[2],
456 hasMother ? kPDecay:kPNoProcess, nt);
7b7936e9 457 // cout << "AliGenUHKM::Generate() pushed on stack with stack index = " << nt
458 // << "; mother index on stack = " << (imo>=0 ? idsOnStack[imo+1] : imo) << endl;
b1c2e580 459
460 idsOnStack[i] = nt;
461 fNprimaries++;
462 KeepTrack(nt);
463
464 origin[0] = origin0[0]+v[0];
465 origin[1] = origin0[1]+v[1];
466 origin[2] = origin0[2]+v[2];
467 imo = nt;
468
7b7936e9 469 trackFlag = kTRUE; // tFlag = kTRUE --> track this one
470 // cout << "AliGenUHKM::Generate() trackFlag = " << trackFlag << endl;
b1c2e580 471
472 PushTrack(trackFlag, imo, kf,
473 p[0], p[1], p[2], energy,
474 origin[0], origin[1], origin[2], iparticle->T(),
475 polar[0], polar[1], polar[2],
476 hasMother ? kPDecay:kPNoProcess, nt);
7b7936e9 477 // cout << "AliGenUHKM::Generate() pushed on stack with stack index = " << nt
478 // << "; mother index on stack = " << imo << endl;
b1c2e580 479 fNprimaries++;
480 KeepTrack(nt);
481 }
482 }
483
484 SetHighWaterMark(fNprimaries);
485
486 TArrayF eventVertex;
487 eventVertex.Set(3);
488 eventVertex[0] = origin0[0];
489 eventVertex[1] = origin0[1];
490 eventVertex[2] = origin0[2];
491
492 // Builds the event header, to be called after each event
493 AliGenEventHeader* header = new AliGenHijingEventHeader("UHKM");
494
495 ((AliGenHijingEventHeader*) header)->SetNProduced(fNprimaries);
496 ((AliGenHijingEventHeader*) header)->SetPrimaryVertex(eventVertex);
497 ((AliGenHijingEventHeader*) header)->SetImpactParameter(0.0);
498 ((AliGenHijingEventHeader*) header)->SetTotalEnergy(0.0);
499 ((AliGenHijingEventHeader*) header)->SetHardScatters(0);
500 ((AliGenHijingEventHeader*) header)->SetParticipants(0, 0);
501 ((AliGenHijingEventHeader*) header)->SetCollisions(0, 0, 0, 0);
502 ((AliGenHijingEventHeader*) header)->SetSpectators(0, 0, 0, 0);
503 ((AliGenHijingEventHeader*) header)->SetReactionPlaneAngle(0);//evrot);
504
505 header->SetPrimaryVertex(fVertex);
506 AddHeader(header);
507 fCollisionGeometry = (AliGenHijingEventHeader*) header;
508
509 delete idsOnStack;
510
511 // gAlice->SetGenEventHeader(header);
512
7b7936e9 513 // printf(" Finish Generate .. %d ..\n",nt);
514 // cout << "AliGenUHKM::Generate() OUT" << endl;
b1c2e580 515}
516
517void AliGenUHKM::Copy(TObject &) const
518{
519 Fatal("Copy","Not implemented!\n");
520}
521
522void AliGenUHKM::SetAllParameters() {
7b7936e9 523 // Forward all input parameters to the TGenerator::TUHKMgen object
524
525 // cout << "AliGenUHKM::SetAllParameters() IN" << endl;
b1c2e580 526
527 fUHKMgen->SetEcms(fHydjetParams.fSqrtS);
528 fUHKMgen->SetBmin(fHydjetParams.fBmin);
529 fUHKMgen->SetBmax(fHydjetParams.fBmax);
530 fUHKMgen->SetAw(fHydjetParams.fAw);
531 fUHKMgen->SetSeed(fHydjetParams.fSeed);
532
533 fUHKMgen->SetChFrzTemperature(fHydjetParams.fT);
534 fUHKMgen->SetMuB(fHydjetParams.fMuB);
535 fUHKMgen->SetMuS(fHydjetParams.fMuS);
536 fUHKMgen->SetMuQ(fHydjetParams.fMuI3);
537 fUHKMgen->SetTauB(fHydjetParams.fTau);
538 fUHKMgen->SetThFrzTemperature(fHydjetParams.fThFO);
539 fUHKMgen->SetMuPionThermal(fHydjetParams.fMu_th_pip);
540
541 fUHKMgen->SetSigmaTau(fHydjetParams.fSigmaTau);
542 fUHKMgen->SetRmaxB(fHydjetParams.fR);
543 fUHKMgen->SetYlMax(fHydjetParams.fYlmax);
544 fUHKMgen->SetEtaRMax(fHydjetParams.fUmax);
545 fUHKMgen->SetMomAsymmPar(fHydjetParams.fDelta);
546 fUHKMgen->SetCoordAsymmPar(fHydjetParams.fEpsilon);
547
548 fUHKMgen->SetGammaS(fHydjetParams.fCorrS);
549 // fUHKMgen->SetHdec(fHydjetParams.fTime);
550 fUHKMgen->SetEtaType(fHydjetParams.fEtaType);
551 fUHKMgen->SetFlagWeakDecay(fHydjetParams.fWeakDecay);
552
553 //PYQUEN parameters
554
555 fUHKMgen->SetPyquenNhsel(fHydjetParams.fNhsel);
556 fUHKMgen->SetPyquenShad(fHydjetParams.fIshad);
557 fUHKMgen->SetPyquenPtmin(fHydjetParams.fPtmin);
558 fUHKMgen->SetPyquenT0(fHydjetParams.fT0);
559 fUHKMgen->SetPyquenTau0(fHydjetParams.fTau0);
560 fUHKMgen->SetPyquenNf(fHydjetParams.fNf);
561 fUHKMgen->SetPyquenIenglu(fHydjetParams.fIenglu);
562 fUHKMgen->SetPyquenIanglu(fHydjetParams.fIanglu);
563
564 fUHKMgen->SetPDGParticleFile(fParticleFilename);
565 fUHKMgen->SetPDGDecayFile(fDecayFilename);
7b7936e9 566 // fUHKMgen->SetUseCharmParticles(fUseCharmParticles);
567 // fUHKMgen->SetMinimumWidth(fMinWidth);
568 // fUHKMgen->SetMaximumWidth(fMaxWidth);
569 // fUHKMgen->SetMinimumMass(fMinMass);
570 // fUHKMgen->SetMaximumMass(fMaxMass);
571
572 cout << "AliGenUHKM::Init() no. stable flagged particles = " << fStableFlagged << endl;
573 for(Int_t i=0; i<fStableFlagged; i++) {
574 cout << "AliGenUHKM::Init() flag no. " << i
575 << " PDG = " << fStableFlagPDG[i]
576 << " flag = " << fStableFlagStatus[i] << endl;
577 fUHKMgen->SetPDGParticleStable(fStableFlagPDG[i], fStableFlagStatus[i]);
578 }
b1c2e580 579
580 cout<<" Print all parameters "<<endl;
581 cout<<" SqrtS = "<<fHydjetParams.fSqrtS<<endl;
582 cout<<" Bmin = "<< fHydjetParams.fBmin<<endl;
583 cout<<" Bmax= "<<fHydjetParams.fBmax<<endl;
584 cout<<" Aw= "<<fHydjetParams.fAw<<endl;
585 cout<<" Seed= "<<fHydjetParams.fSeed<<endl;
586
587 cout<<" ---Stat-model parameters----------- "<<endl;
588
589 cout<<" ChFrzTemperature= "<<fHydjetParams.fT<<endl;
590 cout<<" MuB= "<<fHydjetParams.fMuB<<endl;
591 cout<<" MuS= "<<fHydjetParams.fMuS<<endl;
592 cout<<" MuQ= "<<fHydjetParams.fMuI3<<endl;
593 cout<<" TauB= "<<fHydjetParams.fTau<<endl;
594 cout<<" ThFrzTemperature= "<<fHydjetParams.fThFO<<endl;
595 cout<<" MuPionThermal= "<<fHydjetParams.fMu_th_pip<<endl;
596
597cout<<"-----Volume parameters -------------- "<<endl;
598
599 cout<<" SigmaTau= "<<fHydjetParams.fSigmaTau<<endl;
600 cout<<" RmaxB= "<<fHydjetParams.fR<<endl;
601 cout<<" YlMax= "<<fHydjetParams.fYlmax<<endl;
602 cout<<" EtaRMax= "<<fHydjetParams.fUmax<<endl;
603 cout<<" MomAsymmPar= "<<fHydjetParams.fDelta<<endl;
604 cout<<" CoordAsymmPar= "<<fHydjetParams.fEpsilon<<endl;
605
606cout<<" --------Flags------ "<<endl;
607
608 cout<<" GammaS= "<<fHydjetParams.fCorrS<<endl;
609 // fUHKMgen->SetHdec(fHydjetParams.fTime<<endl;
610 cout<<" EtaType= "<<fHydjetParams.fEtaType<<endl;
611 cout<<" FlagWeakDecay= "<<fHydjetParams.fWeakDecay<<endl;
612
613 cout<<"----PYQUEN parameters---"<<endl;
614
615 cout<<" Nhsel= "<<fHydjetParams.fNhsel<<endl;
616 cout<<" Shad= "<<fHydjetParams.fIshad<<endl;
617 cout<<" Ptmin= "<<fHydjetParams.fPtmin<<endl;
618 cout<<" T0= "<<fHydjetParams.fT0<<endl;
619 cout<<" Tau0= "<<fHydjetParams.fTau0<<endl;
620 cout<<" Nf= "<<fHydjetParams.fNf<<endl;
621 cout<<" Ienglu= "<<fHydjetParams.fIenglu<<endl;
622 cout<<" Ianglu= "<<fHydjetParams.fIanglu<<endl;
623
7b7936e9 624 // cout<<"----PDG table parameters---"<<endl;
b1c2e580 625
7b7936e9 626 // cout<<" UseCharmParticles= "<<fUseCharmParticles<<endl;
627 // cout<<" MinimumWidth= "<<fMinWidth<<endl;
628 // cout<<" MaximumWidth= "<<fMaxWidth<<endl;
629 // cout<<" MinimumMass= "<<fMinMass<<endl;
630 // cout<<" MaximumMass= "<<fMaxMass<<endl;
b1c2e580 631
7b7936e9 632 // cout << "AliGenUHKM::SetAllParameters() OUT" << endl;
b1c2e580 633}
634
635// add the additional PDG codes from UHKM(SHARE table) to ROOT's table
636void AliGenUHKM::CheckPDGTable() {
7b7936e9 637 // Add temporarely all particle definitions from HYDJET++ which miss in the ROOT's PDG tables
638 // to the TDatabasePDG table.
639
640 // cout << "AliGenUHKM::CheckPDGTable() IN" << endl;
b1c2e580 641 //TDabasePDG *rootPDG = TDatabasePDG::Instance(); // ROOT's PDG table
642 DatabasePDG *uhkmPDG = fUHKMgen->PDGInfo(); // UHKM's PDG table
643 TParticlePDG *rootTestParticle;
644 ParticlePDG *uhkmTestParticle;
7b7936e9 645
646 // cout << "particles with good status in UHKM table = " << uhkmPDG->GetNParticles() << endl;
b1c2e580 647 // loop over all particles in the SHARE table
648 for(Int_t i=0; i<uhkmPDG->GetNParticles(); i++) {
7b7936e9 649 // cout << "particle #" << i << " ================" << endl;
b1c2e580 650 // get a particle specie
651 uhkmTestParticle = uhkmPDG->GetPDGParticleByIndex(i);
7b7936e9 652 // cout << "PDG = " << uhkmTestParticle->GetPDG() << endl;
b1c2e580 653 // check if this code exists in ROOT's table
654 rootTestParticle = TDatabasePDG::Instance()->GetParticle(uhkmTestParticle->GetPDG());
655 if(!rootTestParticle) { // if not then add it to the ROOT's PDG database
656
657 TDatabasePDG::Instance()->AddParticle(uhkmTestParticle->GetName(), uhkmTestParticle->GetName(),
658 uhkmTestParticle->GetMass(), uhkmTestParticle->GetElectricCharge(),
659 (uhkmTestParticle->GetWidth()<1e-10 ? kTRUE : kFALSE),
660 uhkmTestParticle->GetWidth(),
661 (Int_t(uhkmTestParticle->GetBaryonNumber())==0 ? "meson" : "baryon"),
662 uhkmTestParticle->GetPDG());
7b7936e9 663 // cout << "Not found in ROOT's PDG database --> added now" << endl;
664 if(uhkmTestParticle->GetWidth()<1e-10)
665 cout << uhkmTestParticle->GetPDG() << "its mass "
666 << TDatabasePDG::Instance()->GetParticle(uhkmTestParticle->GetPDG())->Mass()
667 << TDatabasePDG::Instance()->GetParticle(uhkmTestParticle->GetPDG())->Width() << endl;
b1c2e580 668 }
7b7936e9 669 // else
670 // cout << "Found in ROOT's PDG database --> not added" << endl;
b1c2e580 671 } // end for
7b7936e9 672 // cout << "AliGenUHKM::CheckPDGTable() OUT" << endl;
b1c2e580 673}