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