Coverity fixes
[u/mrichter/AliRoot.git] / TUHKMgen / TUHKMgen.cxx
CommitLineData
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>
31using namespace std;
32
7b7936e9 33//class TObjArray;
34
b1c2e580 35ClassImp(TUHKMgen)
36
37TUHKMgen::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
124 if(strlen(Form("%s/TUHKMgen/UHKM/particles.data", gSystem->Getenv("ALICE_ROOT")))<255)
125 strncpy(fParticleFilename, Form("%s/TUHKMgen/UHKM/particles.data", gSystem->Getenv("ALICE_ROOT")), 256);
126 if(strlen(Form("%s/TUHKMgen/UHKM/tabledecay.txt",gSystem->Getenv("ALICE_ROOT")))<255)
127 strncpy(fDecayFilename, Form("%s/TUHKMgen/UHKM/tabledecay.txt", gSystem->Getenv("ALICE_ROOT")), 256);
b1c2e580 128 for(Int_t i=0; i<500; i++) {
129 fStableFlagPDG[i] = 0;
130 fStableFlagStatus[i] = kFALSE;
131 }
132 fStableFlagged = 0;
7b7936e9 133 // cout << "TUHKMgen::TUHKMgen() OUT" << endl;
b1c2e580 134}
135
136//______________________________________________________________________________
137TUHKMgen::~TUHKMgen()
138{
7b7936e9 139 // destructor, deletes the InitialStateHydjet object
b1c2e580 140 if(fInitialState)
141 delete fInitialState;
b1c2e580 142}
143
144void TUHKMgen::SetAllParametersRHIC()
145{
7b7936e9 146 // Set reasonable input parameters for 0-5% central Au+Au collisions
147 // at 200 GeV at RHIC
b1c2e580 148 SetEcms(200.0); // RHIC top energy
149 SetAw(197); // Au+Au
150 SetIfb(1);
151 SetBfix(0.);
152 SetBmin(0.0); // 0%
153 SetBmax(0.5); // 5%
154 SetChFrzTemperature(0.165); // T_ch = 165 MeV
155 SetMuB(0.0285); // mu_B = 28.5 MeV
156 SetMuS(0.007); // mu_S = 7 MeV
157 SetMuQ(-0.001); // mu_Q = -1 MeV
158 SetThFrzTemperature(0.100); // T_th = 100 MeV
159 SetMuPionThermal(0.053); // mu_th_pion = 53 MeV
160 SetSeed(0); // use UNIX time
161 SetTauB(8.0); // tau = 8 fm/c
162 SetSigmaTau(2.0); // sigma_tau = 2 fm/c
163 SetRmaxB(10.0); // fR = 10 fm
164 SetYlMax(3.3); // fYmax = 3.3
165 SetEtaRMax(1.1); // Umax = 1.1
166 SetMomAsymmPar(0.0); // delta = 0.0
167 SetCoordAsymmPar(0.0); // epsilon = 0.0
168 SetFlagWeakDecay(0.0); // weak decay on (<0 off !!!)
169 SetEtaType(1); // gaus distributed with fYmax dispersion (0 means boost invariant)
170 SetGammaS(1.0); // gammaS = 1.0 (no strangeness canonical suppresion)
171 SetPyquenNhsel(2); // hydro on, jets on, jet quenching on
172 SetPyquenShad(1); // shadowing on (0 off)
173 SetPyquenPtmin(3.4); // ptmin = 3.4 GeV/c
174 SetPyquenT0(0.3); // T0 = 300 MeV
175 SetPyquenTau0(0.4); // tau0 = 0.4 fm/c
176 SetPyquenNf(2); // 2 flavours
177 SetPyquenIenglu(0); // radiative and collisional energy loss
178 SetPyquenIanglu(0); // small gluon angular distribution
179}
180
181void TUHKMgen::SetAllParametersLHC()
182{
7b7936e9 183 // Set reasonable input parameters for 0-5% Pb+Pb collisions at 5.5 TeV at LHC
b1c2e580 184 SetEcms(5500.0); // LHC
185 SetAw(207); // Pb+Pb
186 SetIfb(1);
187 SetBfix(0.); // 0
188 SetBmin(0.0); // 0%
189 SetBmax(0.5); // 5%
190 SetChFrzTemperature(0.170); // T_ch = 170 MeV
191 SetMuB(0.0); // mu_B = 0 MeV
192 SetMuS(0.0); // mu_S = 0 MeV
193 SetMuQ(0.0); // mu_Q = 0 MeV
194 SetThFrzTemperature(0.130); // T_th = 130 MeV
195 SetMuPionThermal(0.0); // mu_th_pion = 0 MeV
196 SetSeed(0); // use UNIX time
197 SetTauB(10.0); // tau = 10 fm/c
198 SetSigmaTau(3.0); // sigma_tau = 3 fm/c
199 SetRmaxB(11.0); // fR = 11 fm
200 SetYlMax(4.0); // fYmax = 4.0
201 SetEtaRMax(1.1); // Umax = 1.1
202 SetMomAsymmPar(0.0); // delta = 0.0
203 SetCoordAsymmPar(0.0); // epsilon = 0.0
204 SetFlagWeakDecay(0); // weak decay on (<0 off !!!)
205 SetEtaType(1); // gaus distributed with fYmax dispersion (0 means boost invariant)
206 SetGammaS(1.0); // gammaS = 1.0 (no strangeness canonical suppresion)
207 SetPyquenNhsel(2); // hydro on, jets on, jet quenching on
208 SetPyquenShad(1); // shadowing on (0 off)
209 SetPyquenPtmin(7.0); // ptmin = 7.0 GeV/c
210 SetPyquenT0(0.8); // T0 = 800 MeV
211 SetPyquenTau0(0.1); // tau0 = 0.4 fm/c
212 SetPyquenNf(0); // 0 flavours
213 SetPyquenIenglu(0); // radiative and collisional energy loss
214 SetPyquenIanglu(0); // small gluon angular distribution
215}
216
7b7936e9 217TObjArray* TUHKMgen::ImportParticles(const Option_t *)
b1c2e580 218{
7b7936e9 219 // Function overloading the TGenerator::ImportParticles() member function.
220 // The particles from the local particle list (fSecondariesList) are
221 // forwarded to the TGenerator::fParticles
b1c2e580 222 fParticles->Clear();
b1c2e580 223 Int_t nump = 0;
224 LPIT_t it,e;
b1c2e580 225
3fa37a65 226 for(it = fSecondariesList.begin(), e = fSecondariesList.end(); it != e; it++) {
b1c2e580 227 TVector3 pos(it->Pos().Vect());
228 TVector3 mom(it->Mom().Vect());
229 Float_t m1 = it->TableMass();
230 Int_t im1 = it->GetMother();
231 Int_t im2 = -1;
232 Int_t id1 = -1;
233 Int_t id2 = -1;
b1c2e580 234
b1c2e580 235 Int_t type = it->GetType(); // 0-hydro, 1-jets
236
237 if (im1> -1) {
3fa37a65 238 TParticle *mother = (TParticle*) (fParticles->UncheckedAt(im1+1));
b1c2e580 239 mother->SetLastDaughter(nump);
240 if (mother->GetFirstDaughter()==-1)
3fa37a65 241 mother->SetFirstDaughter(nump+1);
b1c2e580 242 }
243
244 nump++;
245
246 TParticle* p = new TParticle(it->Encoding(), type, //pdg,stat
247 im1, im2, id1, id2, //m1,m2,d1,d2
248 mom[0], mom[1], mom[2], TMath::Sqrt(mom.Mag2() + m1 * m1), //px,py,pz,e
7b7936e9 249 pos[0]*1.e-13, pos[1]*1.e-13, pos[2]*1.e-13,
250 it->T()*1.e-13/3e10); //x,y,z,t
b1c2e580 251
252 p->SetUniqueID(nump);
253 fParticles->Add(p);
254 } //end for
255
b1c2e580 256 fAllocator.FreeList(fSecondariesList);
3fa37a65 257
b1c2e580 258 return fParticles;
259}
260
7b7936e9 261Int_t TUHKMgen::ImportParticles(TClonesArray *particles, const Option_t* option)
b1c2e580 262{
7b7936e9 263 // Function overloading the TGenerator::ImportParticles() member function.
264 // The particles from the local particle list (fSecondariesList) are
265 // forwarded to the TGenerator::fParticles
25c1936c 266 option = option; // just to avoid the warning
3fa37a65 267
268
b1c2e580 269 if(particles==0) return 0;
270 TClonesArray &particlesR=*particles;
271 particlesR.Clear();
b1c2e580 272
273 Int_t numprim,numsec; numprim=numsec=0;
274 Int_t nump = 0;
275 LPIT_t it,e;
7b7936e9 276
3fa37a65 277 for(it = fSecondariesList.begin(), e = fSecondariesList.end(); it != e; it++) {
278 TVector3 pos(it->Pos().Vect());
279 TVector3 mom(it->Mom().Vect());
280 Float_t m1 = it->TableMass();
281 Int_t im1 = it->GetMother();
b1c2e580 282 Int_t im2 = -1;
283 Int_t id1 = -1;
284 Int_t id2 = -1;
b1c2e580 285
3fa37a65 286 Int_t type = it->GetType();
287
b1c2e580 288 if (im1> -1) {
7fecd07a 289 // particle not a primary -> set the daughter indexes for the mother particle"<< endl;
3fa37a65 290 TParticle *mother = (TParticle*) (particlesR.UncheckedAt(im1));
b1c2e580 291 mother->SetLastDaughter(nump);
3fa37a65 292 if(mother->GetFirstDaughter()==-1)
293 mother->SetFirstDaughter(nump);
7fecd07a 294 } else
295 ++numprim;
3fa37a65 296
b1c2e580 297 new (particlesR[nump]) TParticle(it->Encoding(), type, //pdg,stat
298 im1, im2, id1, id2, //m1,m2,d1,d2
299 mom[0], mom[1], mom[2], TMath::Sqrt(mom.Mag2() + m1 * m1), //px,py,pz,e
7b7936e9 300 pos[0]*1.e-13, pos[1]*1.e-13, pos[2]*1.e-13,
301 it->T()*1.e-13/3e10); //x,y,z,t
302
b1c2e580 303 particlesR[nump]->SetUniqueID(nump);
3fa37a65 304 nump++;
b1c2e580 305 numsec++;
306 }//end for
307
3fa37a65 308 fSecondariesList.clear();
309 printf("Scan and add prim %d sec %d and all %d particles\n",
310 numprim,numsec,nump);
b1c2e580 311 return nump;
312}
313
314//______________________________________________________________________________
315void TUHKMgen::Initialize()
316{
7b7936e9 317 // Function overloading the TGenerator::Initialize() member function.
318 // The Monte-Carlo model is initialized (input parameters are transmitted,
319 // particle list and decay channels are loaded, average multiplicities are calculated, etc.)
3fa37a65 320
321 fInitialState = new InitialStateHydjet();
322 SetAllParameters();
323 fInitialState->LoadPDGInfo();
b1c2e580 324 // set the stable flags
3fa37a65 325 for(Int_t i=0; i<fStableFlagged; i++)
326 fInitialState->SetPDGParticleStable(fStableFlagPDG[i], fStableFlagStatus[i]);
b1c2e580 327
3fa37a65 328 if(!fInitialState->MultIni())
7b7936e9 329 Error("TUHKMgen::Initialize()", "Bad status return from MultIni(). Check it out!! \n"); //
b1c2e580 330}
331
786056a2 332void TUHKMgen::Print(const Option_t*) const
b1c2e580 333{
b7dbe0fd 334 cout << "**********************************************************************************" << endl;
335 cout << "* UHKM Generator interface to ROOT::TGenerator *" << endl;
336 cout << "* Documentation: *" << endl;
337 cout << "* I.P.Lokhtin, L.V.Malinina, S.V.Petrushanko, A.M.Snigirev, I.Arsene, K.Tywoniuk *" << endl;
338 cout << "* Comput.Phys.Commun.180:779-799, 2009 *" << endl;
339 cout << "**********************************************************************************" << endl;
b1c2e580 340}
341
342void TUHKMgen::GenerateEvent()
343{
7b7936e9 344 // Member function overloading the TGenerator::GenerateEvent()
345 // The HYDJET++ model is run and the particle lists (fSourceList and fSecondariesList) are filled
7b7936e9 346
3fa37a65 347 fInitialState->Initialize(fSecondariesList, fAllocator);
b1c2e580 348
3fa37a65 349 if(fSecondariesList.empty())Error("TUHKMgen::GenerateEvent()", "Source particle list empty after fireball initialization!! \n"); //
b1c2e580 350
351 // Run the decays
7b7936e9 352
3fa37a65 353 if(fInitialState->RunDecays())
354 fInitialState->Evolve(fSecondariesList, fAllocator, fInitialState->GetWeakDecayLimit());
355
b1c2e580 356}
357
358void TUHKMgen::SetAllParameters() {
7b7936e9 359 // forward all model input parameters to the InitialStateHydjet object
360 // which will handle all the Monte-Carlo simulation using HYDJET++ model
361
b1c2e580 362 fInitialState->fParams.fSqrtS = fHydjetParams.fSqrtS;
363 fInitialState->fParams.fAw = fHydjetParams.fAw;
364 fInitialState->fParams.fIfb = fHydjetParams.fIfb ;
365 fInitialState->fParams.fBfix = fHydjetParams.fBfix;
366 fInitialState->fParams.fBmin = fHydjetParams.fBmin;
367 fInitialState->fParams.fBmax = fHydjetParams.fBmax;
368 fInitialState->fParams.fSeed = fHydjetParams.fSeed;
369 fInitialState->fParams.fT = fHydjetParams.fT;
370 fInitialState->fParams.fMuB = fHydjetParams.fMuB;
371 fInitialState->fParams.fMuS = fHydjetParams.fMuS;
372 fInitialState->fParams.fMuI3 = fHydjetParams.fMuI3;
373 fInitialState->fParams.fThFO = fHydjetParams.fThFO;
374 fInitialState->fParams.fMu_th_pip = fHydjetParams.fMu_th_pip;
375
376 fInitialState->fParams.fTau = fHydjetParams.fTau;
377 fInitialState->fParams.fSigmaTau = fHydjetParams.fSigmaTau;
378 fInitialState->fParams.fR = fHydjetParams.fR;
379 fInitialState->fParams.fYlmax = fHydjetParams.fYlmax;
380 fInitialState->fParams.fUmax = fHydjetParams.fUmax;
381 fInitialState->fParams.fDelta = fHydjetParams.fDelta;
382 fInitialState->fParams.fEpsilon = fHydjetParams.fEpsilon;
383
b1c2e580 384 fInitialState->fParams.fWeakDecay = fHydjetParams.fWeakDecay;
3fa37a65 385 fInitialState->fParams.fDecay = 1; // run decays
b1c2e580 386 fInitialState->fParams.fCorrS = fHydjetParams.fCorrS;
387 fInitialState->fParams.fEtaType = fHydjetParams.fEtaType;
388
389 fInitialState->fParams.fPtmin = fHydjetParams.fPtmin;
390 fInitialState->fParams.fNhsel = fHydjetParams.fNhsel;
391 fInitialState->fParams.fIshad = fHydjetParams.fIshad;
392 fInitialState->fParams.fT0 = fHydjetParams.fT0;
393 fInitialState->fParams.fTau0 = fHydjetParams.fTau0;
394 fInitialState->fParams.fNf = fHydjetParams.fNf;
395 fInitialState->fParams.fIenglu = fHydjetParams.fIenglu;
396 fInitialState->fParams.fIanglu = fHydjetParams.fIanglu;
397
398 fInitialState->SetPDGParticleFilename(fParticleFilename);
399 fInitialState->SetPDGDecayFilename(fDecayFilename);
7b7936e9 400 // fInitialState->SetUseCharmParticles(fUseCharmParticles);
401 // fInitialState->SetWidthRange(fMinWidth, fMaxWidth);
402 // fInitialState->SetMassRange(fMinMass, fMaxMass);
403 // for(Int_t i=0; i<fStableFlagged; i++)
404 // fInitialState->SetPDGParticleStable(fStableFlagPDG[i], fStableFlagStatus[i]);
405 // cout << "TUHKMgen::SetAllParameters() OUT" << endl;
b1c2e580 406}