71c1902d221d96f9b4f88664a5f40183aad03ad2
[u/mrichter/AliRoot.git] / TUHKMgen / AliGenUHKM.cxx
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
17 #include <iostream>
18 #include <string>
19
20 #include "TUHKMgen.h"
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>
32 #include <TDatabasePDG.h>
33 #include <TSystem.h>
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
43 using namespace std;
44
45
46 ClassImp(AliGenUHKM)
47
48 //_______________________________________
49 AliGenUHKM::AliGenUHKM()
50   :AliGenMC(),
51    fTrials(0),
52    fUHKMgen(0),
53    fHydjetParams(),
54    fStableFlagged(0)
55 {
56   // Default constructor setting up default reasonable parameter values
57   // for central Pb+Pb collisions at 5.5TeV 
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 */
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);
126   for(Int_t i=0; i<500; i++) {
127     fStableFlagPDG[i] = 0;
128     fStableFlagStatus[i] = kFALSE;
129   }
130   fStableFlagged = 0;
131
132 }
133
134 //_______________________________________
135 AliGenUHKM::AliGenUHKM(Int_t npart)
136   :AliGenMC(npart),
137    fTrials(0),
138    fUHKMgen(0),
139    fHydjetParams(),
140    fStableFlagged(0)
141 {
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
145
146   fName = "UHKM";
147   fTitle= "Particle Generator using UHKM 3.0";
148
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
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);
218   for(Int_t i=0; i<500; i++) {
219     fStableFlagPDG[i] = 0;
220     fStableFlagStatus[i] = kFALSE;
221   }
222   fStableFlagged = 0;  
223
224 }
225
226 //__________________________________________
227 AliGenUHKM::~AliGenUHKM()
228 {
229   // Destructor, do nothing
230   //  delete fParticles;
231 }
232
233 void AliGenUHKM::SetAllParametersRHIC() 
234 {
235   // Set reasonable default parameters for 0-5% central Au+Au collisions
236   // at 200 GeV at RHIC
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
255   //  SetFlagWeakDecay(0);           // weak decay on (<0 off !!!)
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
268 void AliGenUHKM::SetAllParametersLHC()
269 {
270   // Set reasonable default parameters for 0-5% central Pb+Pb collisions
271   // at 5.5 TeV at LHC
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
290   //  SetFlagWeakDecay(0);           // weak decay on (<0 off !!!)
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 //_________________________________________
304 void AliGenUHKM::Init()
305 {
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.)
311
312   SetMC(new TUHKMgen());
313   fUHKMgen = (TUHKMgen*) fMCEvGen;
314   SetAllParameters();
315   
316   AliGenMC::Init();
317   fUHKMgen->Initialize();
318   CheckPDGTable();    
319   fUHKMgen->Print();  
320 }
321
322 //________________________________________
323 void AliGenUHKM::Generate()
324 {
325   // Generate one HYDJET++ event, get the output and push particles further 
326   // to AliRoot's stack 
327
328   Float_t polar[3] = {0,0,0};
329   Float_t origin[3]   = {0,0,0};
330   Float_t origin0[3]  = {0,0,0};
331   Float_t p[3];
332   Float_t v[3];
333   Float_t mass=0.0, energy=0.0;
334
335   Vertex();
336   for(Int_t j=0; j<3; j++) origin0[j] = fVertex[j];
337
338   // Generate the event and import particles
339   fUHKMgen->GenerateEvent();
340   fUHKMgen->ImportParticles(&fParticles,"All");
341   Int_t np = fParticles.GetEntriesFast();
342   Int_t nt  =  0;
343
344
345   // Handle the IDs of particles on the stack  
346   Int_t* idsOnStack = new Int_t[np];
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   }
352   
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
361   for(Int_t i=0; i<np; i++) {
362     iparticle = (TParticle*)fParticles.At(i);
363     Int_t kf = iparticle->GetPdgCode();
364     Bool_t hasMother = (iparticle->GetFirstMother() >= 0);
365     Bool_t hasDaughter = (iparticle->GetNDaughters() > 0);
366     
367     if(hasDaughter) {
368       // This particle has decayed
369       // It will not be tracked
370       // Add it only once with coordinates not
371       // smeared with primary vertex position
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       
380       mass = TDatabasePDG::Instance()->GetParticle(kf)->Mass();
381       energy = sqrt(mass*mass + p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
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);
388       v[2] = iparticle->Vz();
389
390       Float_t time = iparticle->T();
391       
392       Int_t imo = -1;
393       if(hasMother) {
394         imo = iparticle->GetFirstMother(); //index of mother particle in fParticles
395       } // if has mother
396       Bool_t trackFlag = kFALSE;   // tFlag is kFALSE --> do not track the particle
397       
398       PushTrack(trackFlag, (imo>=0 ? idsOnStack[imo] : imo), kf,
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);
403       idsOnStack[i] = nt;
404       
405       fNprimaries++;
406       KeepTrack(nt);
407     }
408     else {
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
414       // Second time with event-wide c0ordinates and vertex smearing
415       //   this one will be tracked
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
424       energy = sqrt(mass*mass + p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
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();      
432       
433       Int_t type    = iparticle->GetStatusCode(); // 1-from jet / 0-from hydro 
434       Int_t coeffT=1;
435       if(type==1) coeffT=-1; //to separate particles from jets
436       
437       
438       Int_t imo = -1;
439       
440       if(hasMother) {
441         imo = iparticle->GetFirstMother();
442       } // if has mother
443       
444       Bool_t trackFlag = kFALSE;  // tFlag = kFALSE --> do not track this one, its for femtoscopy
445       PushTrack(trackFlag, (imo>=0 ? idsOnStack[imo] : imo), kf,
446                 p[0], p[1], p[2], energy,
447                 v[0], v[1], v[2], (iparticle->T())*coeffT,    // freeze-out time is negative if the particle comes from jet
448                 polar[0], polar[1], polar[2],
449                 hasMother ? kPDecay:kPNoProcess, nt);
450       
451       idsOnStack[i] = nt;
452       fNprimaries++;
453       KeepTrack(nt);
454       
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       
460       trackFlag = fTrackIt;    // Track this particle, unless otherwise specified by fTrackIt
461       
462       PushTrack(trackFlag, imo, kf,
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);
467       
468       fNprimaries++;
469       KeepTrack(nt);
470     }
471   }
472   
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");
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
490
491   ((AliGenHijingEventHeader*) header)->SetNProduced(fNprimaries);
492   ((AliGenHijingEventHeader*) header)->SetPrimaryVertex(eventVertex);
493   ((AliGenHijingEventHeader*) header)->SetImpactParameter(b);
494   ((AliGenHijingEventHeader*) header)->SetTotalEnergy(0.0);
495   ((AliGenHijingEventHeader*) header)->SetHardScatters(0);
496   ((AliGenHijingEventHeader*) header)->SetParticipants(Int_t(npart), 0);
497   ((AliGenHijingEventHeader*) header)->SetCollisions(Int_t(nbin), 0, 0, 0);
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
505   delete [] idsOnStack;
506
507 }
508
509 void AliGenUHKM::Copy(TObject &) const
510 {
511   Fatal("Copy","Not implemented!\n");
512 }
513
514 void AliGenUHKM::SetAllParameters() {
515   // Forward all input parameters to the TGenerator::TUHKMgen object
516
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);
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);
555   //  fUHKMgen->SetUseCharmParticles(fUseCharmParticles);
556   //  fUHKMgen->SetMinimumWidth(fMinWidth);
557   //  fUHKMgen->SetMaximumWidth(fMaxWidth);
558   //  fUHKMgen->SetMinimumMass(fMinMass);
559   //  fUHKMgen->SetMaximumMass(fMaxMass);
560
561   //  cout << "AliGenUHKM::Init() no. stable flagged particles = " << fStableFlagged << endl;
562   for(Int_t i=0; i<fStableFlagged; i++) {
563     //    cout << "AliGenUHKM::Init() flag no. " << i
564     //         << " PDG = " << fStableFlagPDG[i]
565     //         << " flag = " << fStableFlagStatus[i] << endl;
566     fUHKMgen->SetPDGParticleStable(fStableFlagPDG[i], fStableFlagStatus[i]);
567   }
568
569    cout<<" Print all parameters "<<endl;
570    cout<<" SqrtS = "<<fHydjetParams.fSqrtS<<endl;
571    cout<<" Bmin  = "<< fHydjetParams.fBmin<<endl;
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
586 cout<<"-----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
595 cout<<" --------Flags------ "<<endl;
596
597  cout<<" GammaS= "<<fHydjetParams.fCorrS<<endl;
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
612   //  cout<<"----PDG table parameters---"<<endl;
613   
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;
619
620   //  cout << "AliGenUHKM::SetAllParameters() OUT" << endl;
621 }
622
623 // add the additional PDG codes from UHKM(SHARE table) to ROOT's table
624 void AliGenUHKM::CheckPDGTable() {
625   // Add temporarely all particle definitions from HYDJET++ which miss in the ROOT's PDG tables
626   // to the TDatabasePDG table.
627
628   DatabasePDG *uhkmPDG = fUHKMgen->PDGInfo();              // UHKM's PDG table
629   TParticlePDG *rootTestParticle;
630   ParticlePDG *uhkmTestParticle;
631   
632   // loop over all particles in the SHARE table
633   for(Int_t i=0; i<uhkmPDG->GetNParticles(); i++) {
634     // get a particle specie
635     uhkmTestParticle = uhkmPDG->GetPDGParticleByIndex(i);
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());    
646       if(uhkmTestParticle->GetWidth()<1e-10) 
647         cout << uhkmTestParticle->GetPDG() << " with mass "
648              << TDatabasePDG::Instance()->GetParticle(uhkmTestParticle->GetPDG())->Mass()
649              << TDatabasePDG::Instance()->GetParticle(uhkmTestParticle->GetPDG())->Width() << endl;  
650     }
651   }  // end for
652 }