Coding conventions (Ionut)
[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*/
123
124 strcpy(fParticleFilename, Form("%s/TUHKMgen/UHKM/particles.data", gSystem->Getenv("ALICE_ROOT")));
125 strcpy(fDecayFilename, Form("%s/TUHKMgen/UHKM/tabledecay.txt", gSystem->Getenv("ALICE_ROOT")));
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//______________________________________________________________________________
135TUHKMgen::~TUHKMgen()
136{
7b7936e9 137 // destructor, deletes the InitialStateHydjet object
b1c2e580 138 if(fInitialState)
139 delete fInitialState;
b1c2e580 140}
141
142void 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
179void 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 215TObjArray* 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 259Int_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) {
3fa37a65 287 // particle not a primary -> set the daughter indexes for the mother particle"<< endl;
288 TParticle *mother = (TParticle*) (particlesR.UncheckedAt(im1));
b1c2e580 289 mother->SetLastDaughter(nump);
3fa37a65 290 if(mother->GetFirstDaughter()==-1)
291 mother->SetFirstDaughter(nump);
b1c2e580 292 }
293
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//______________________________________________________________________________
313void 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 330void TUHKMgen::Print(const Option_t*) const
b1c2e580 331{
332 cout << "TUHKMgen::Print() method not implemented yet!!" << endl;
333}
334
335void TUHKMgen::GenerateEvent()
336{
7b7936e9 337 // Member function overloading the TGenerator::GenerateEvent()
338 // The HYDJET++ model is run and the particle lists (fSourceList and fSecondariesList) are filled
7b7936e9 339
3fa37a65 340 fInitialState->Initialize(fSecondariesList, fAllocator);
b1c2e580 341
3fa37a65 342 if(fSecondariesList.empty())Error("TUHKMgen::GenerateEvent()", "Source particle list empty after fireball initialization!! \n"); //
b1c2e580 343
344 // Run the decays
7b7936e9 345
3fa37a65 346 if(fInitialState->RunDecays())
347 fInitialState->Evolve(fSecondariesList, fAllocator, fInitialState->GetWeakDecayLimit());
348
b1c2e580 349}
350
351void TUHKMgen::SetAllParameters() {
7b7936e9 352 // forward all model input parameters to the InitialStateHydjet object
353 // which will handle all the Monte-Carlo simulation using HYDJET++ model
354
b1c2e580 355 fInitialState->fParams.fSqrtS = fHydjetParams.fSqrtS;
356 fInitialState->fParams.fAw = fHydjetParams.fAw;
357 fInitialState->fParams.fIfb = fHydjetParams.fIfb ;
358 fInitialState->fParams.fBfix = fHydjetParams.fBfix;
359 fInitialState->fParams.fBmin = fHydjetParams.fBmin;
360 fInitialState->fParams.fBmax = fHydjetParams.fBmax;
361 fInitialState->fParams.fSeed = fHydjetParams.fSeed;
362 fInitialState->fParams.fT = fHydjetParams.fT;
363 fInitialState->fParams.fMuB = fHydjetParams.fMuB;
364 fInitialState->fParams.fMuS = fHydjetParams.fMuS;
365 fInitialState->fParams.fMuI3 = fHydjetParams.fMuI3;
366 fInitialState->fParams.fThFO = fHydjetParams.fThFO;
367 fInitialState->fParams.fMu_th_pip = fHydjetParams.fMu_th_pip;
368
369 fInitialState->fParams.fTau = fHydjetParams.fTau;
370 fInitialState->fParams.fSigmaTau = fHydjetParams.fSigmaTau;
371 fInitialState->fParams.fR = fHydjetParams.fR;
372 fInitialState->fParams.fYlmax = fHydjetParams.fYlmax;
373 fInitialState->fParams.fUmax = fHydjetParams.fUmax;
374 fInitialState->fParams.fDelta = fHydjetParams.fDelta;
375 fInitialState->fParams.fEpsilon = fHydjetParams.fEpsilon;
376
b1c2e580 377 fInitialState->fParams.fWeakDecay = fHydjetParams.fWeakDecay;
3fa37a65 378 fInitialState->fParams.fDecay = 1; // run decays
b1c2e580 379 fInitialState->fParams.fCorrS = fHydjetParams.fCorrS;
380 fInitialState->fParams.fEtaType = fHydjetParams.fEtaType;
381
382 fInitialState->fParams.fPtmin = fHydjetParams.fPtmin;
383 fInitialState->fParams.fNhsel = fHydjetParams.fNhsel;
384 fInitialState->fParams.fIshad = fHydjetParams.fIshad;
385 fInitialState->fParams.fT0 = fHydjetParams.fT0;
386 fInitialState->fParams.fTau0 = fHydjetParams.fTau0;
387 fInitialState->fParams.fNf = fHydjetParams.fNf;
388 fInitialState->fParams.fIenglu = fHydjetParams.fIenglu;
389 fInitialState->fParams.fIanglu = fHydjetParams.fIanglu;
390
391 fInitialState->SetPDGParticleFilename(fParticleFilename);
392 fInitialState->SetPDGDecayFilename(fDecayFilename);
7b7936e9 393 // fInitialState->SetUseCharmParticles(fUseCharmParticles);
394 // fInitialState->SetWidthRange(fMinWidth, fMaxWidth);
395 // fInitialState->SetMassRange(fMinMass, fMaxMass);
396 // for(Int_t i=0; i<fStableFlagged; i++)
397 // fInitialState->SetPDGParticleStable(fStableFlagPDG[i], fStableFlagStatus[i]);
398 // cout << "TUHKMgen::SetAllParameters() OUT" << endl;
b1c2e580 399}