Removing extra symbols
[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 <TUHKMgen.h>
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>
42 using namespace std;
43
44 ClassImp(AliGenUHKM)
45
46 //_______________________________________
47 AliGenUHKM::AliGenUHKM()
48   :AliGenMC(),
49    fTrials(0),
50    fUHKMgen(0),
51    fUseCharmParticles(kFALSE),
52    fMinWidth(0.0),
53    fMaxWidth(10.0),
54    fMinMass(0.0001),
55    fMaxMass(10.0)
56 {
57   cout << "AliGenUHKM::AliGenUHKM() IN" << endl;
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   strcpy(fParticleFilename, Form("%s/TUHKMgen/UHKM/particles.data", gSystem->Getenv("ALICE_ROOT")));
123   strcpy(fDecayFilename, Form("%s/TUHKMgen/UHKM/tabledecay.txt", gSystem->Getenv("ALICE_ROOT")));
124   for(Int_t i=0; i<500; i++) {
125     fStableFlagPDG[i] = 0;
126     fStableFlagStatus[i] = kFALSE;
127   }
128   fStableFlagged = 0;
129
130   cout << "AliGenUHKM::AliGenUHKM() OUT" << endl;
131 }
132
133 //_______________________________________
134 AliGenUHKM::AliGenUHKM(Int_t npart)
135   :AliGenMC(npart),
136    fTrials(0),
137    fUHKMgen(0),
138    fUseCharmParticles(kFALSE),
139    fMinWidth(0.0),
140    fMaxWidth(10.0),
141    fMinMass(0.0001),
142    fMaxMass(10.0)
143 {
144   cout << "AliGenUHKM::AliGenUHKM(Int_t) IN" << endl;
145   fName = "UHKM";
146   fTitle= "Particle Generator using UHKM 3.0";
147
148   // Constructor specifying the size of the particle table
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   strcpy(fParticleFilename, Form("%s/TUHKMgen/UHKM/particles.data", gSystem->Getenv("ALICE_ROOT")));
215   strcpy(fDecayFilename, Form("%s/TUHKMgen/UHKM/tabledecay.txt", gSystem->Getenv("ALICE_ROOT")));
216   for(Int_t i=0; i<500; i++) {
217     fStableFlagPDG[i] = 0;
218     fStableFlagStatus[i] = kFALSE;
219   }
220   fStableFlagged = 0;  
221
222   cout << "AliGenUHKM::AliGenUHKM(Int_t) OUT" << endl;
223 }
224
225 //__________________________________________
226 AliGenUHKM::~AliGenUHKM()
227 {}
228
229 void AliGenUHKM::SetAllParametersRHIC() 
230 {
231   SetEcms(200.0);                // RHIC top energy
232   SetAw(197);                    // Au+Au
233   SetBmin(0.0);                  // 0%
234   SetBmax(0.5);                  // 5%
235   SetChFrzTemperature(0.165);    // T_ch = 165 MeV
236   SetMuB(0.0285);                // mu_B = 28.5 MeV
237   SetMuS(0.007);                 // mu_S = 7 MeV
238   SetMuQ(-0.001);                // mu_Q = -1 MeV
239   SetThFrzTemperature(0.100);    // T_th = 100 MeV
240   SetMuPionThermal(0.053);       // mu_th_pion = 53 MeV 
241   SetSeed(0);                    // use UNIX time
242   SetTauB(8.0);                  // tau = 8 fm/c
243   SetSigmaTau(2.0);              // sigma_tau = 2 fm/c
244   SetRmaxB(10.0);                // fR = 10 fm
245   SetYlMax(3.3);                 // fYmax = 3.3
246   SetEtaRMax(1.1);               // Umax = 1.1
247   SetMomAsymmPar(0.0);           // delta = 0.0
248   SetCoordAsymmPar(0.0);         // epsilon = 0.0
249   SetFlagWeakDecay(0);           // weak decay on (<0 off !!!)
250   SetEtaType(1);                 // gaus distributed with fYmax dispersion (0 means boost invariant)
251   SetGammaS(1.0);                // gammaS = 1.0 (no strangeness canonical suppresion)
252   SetPyquenNhsel(2);             // hydro on, jets on, jet quenching on
253   SetPyquenShad(1);              // shadowing on (0 off)
254   SetPyquenPtmin(3.4);           // ptmin = 3.4 GeV/c
255   SetPyquenT0(0.3);              // T0 = 300 MeV
256   SetPyquenTau0(0.4);            // tau0 = 0.4 fm/c
257   SetPyquenNf(2);                // 2 flavours
258   SetPyquenIenglu(0);            // radiative and collisional energy loss
259   SetPyquenIanglu(0);            // small gluon angular distribution
260 }
261
262 void AliGenUHKM::SetAllParametersLHC()
263 {
264   SetEcms(5500.0);               // LHC
265   SetAw(207);                    // Pb+Pb
266   SetBmin(0.0);                  // 0%
267   SetBmax(0.5);                  // 5%
268   SetChFrzTemperature(0.170);    // T_ch = 170 MeV
269   SetMuB(0.0);                   // mu_B = 0 MeV
270   SetMuS(0.0);                   // mu_S = 0 MeV
271   SetMuQ(0.0);                   // mu_Q = 0 MeV
272   SetThFrzTemperature(0.130);    // T_th = 130 MeV
273   SetMuPionThermal(0.0);         // mu_th_pion = 0 MeV
274   SetSeed(0);                    // use UNIX time
275   SetTauB(10.0);                 // tau = 10 fm/c
276   SetSigmaTau(3.0);              // sigma_tau = 3 fm/c
277   SetRmaxB(11.0);                // fR = 11 fm
278   SetYlMax(4.0);                 // fYmax = 4.0
279   SetEtaRMax(1.1);               // Umax = 1.1
280   SetMomAsymmPar(0.0);           // delta = 0.0
281   SetCoordAsymmPar(0.0);         // epsilon = 0.0
282   SetFlagWeakDecay(0);           // weak decay on (<0 off !!!)
283   SetEtaType(1);                 // gaus distributed with fYmax dispersion (0 means boost invariant)
284   SetGammaS(1.0);                // gammaS = 1.0 (no strangeness canonical suppresion)
285   SetPyquenNhsel(2);             // hydro on, jets on, jet quenching on
286   SetPyquenShad(1);              // shadowing on (0 off)
287   SetPyquenPtmin(7.0);           // ptmin = 7.0 GeV/c
288   SetPyquenT0(0.8);              // T0 = 800 MeV
289   SetPyquenTau0(0.1);            // tau0 = 0.4 fm/c
290   SetPyquenNf(0);                // 0 flavours
291   SetPyquenIenglu(0);            // radiative and collisional energy loss
292   SetPyquenIanglu(0);            // small gluon angular distribution
293 }
294
295 //_________________________________________
296 void AliGenUHKM::Init()
297 {
298   cout << "AliGenUHKM::Init() IN" << endl;
299
300   SetMC(new TUHKMgen());
301   fUHKMgen = (TUHKMgen*) fMCEvGen;
302   SetAllParameters();
303
304   AliGenMC::Init();
305
306   fUHKMgen->Initialize();
307   CheckPDGTable();
308
309   fUHKMgen->Print();
310   cout << "AliGenUHKM::Init() OUT" << endl;
311 }
312
313
314
315 //________________________________________
316 void AliGenUHKM::Generate()
317 {
318   cout << "AliGenUHKM::Generate() IN" << endl;
319   Float_t polar[3] = {0,0,0};
320   Float_t origin[3]   = {0,0,0};
321   Float_t origin0[3]  = {0,0,0};
322   Float_t v[3];
323   Float_t mass, energy;
324
325   Vertex();
326   for(Int_t j=0; j<3; j++) origin0[j] = fVertex[j];
327
328   fTrials = 0;
329  // cout << "AliGenUHKM::Generate() fTrials = " << fTrials << endl;
330
331   Int_t nt  = 0;
332
333   fUHKMgen->GenerateEvent();
334   fTrials++;
335
336   fUHKMgen->ImportParticles(&fParticles,"All");
337
338   Int_t np = fParticles.GetEntriesFast();
339   cout << "AliGenUHKM::Generate() GetEntries  " <<np<< endl;
340
341
342   Int_t* idsOnStack = new Int_t[np];
343   Int_t* newPos     = new Int_t[np];
344   for(Int_t i=0; i<np; i++) newPos[i] = i;
345
346   //_________ Loop for particle selection
347   //  for(Int_t i=1; i<np; i++) {
348   for(Int_t i=1; i<np; i++) {
349     // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
350     // the particle indexes are 0 based but fParticles is a 1 based array
351     // -1 is the trivial code (when it does not exist)
352     // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
353     TParticle *iparticle = (TParticle*)fParticles.At(i);
354 //    cout << "AliGenUHKM::Generate() particle #" << i << " in fParticles *********************"<< endl;
355
356     Int_t kf = iparticle->GetPdgCode();
357 //    cout << "AliGenUHKM::Generate() PDG = " << kf << endl;
358
359     Bool_t hasMother = (iparticle->GetFirstMother() >= 0);
360
361 //    cout << "AliGenUHKM::Generate() mother index in fParticles = " 
362 //       << (iparticle->GetFirstMother()==-1 ? -1 : iparticle->GetFirstMother()+1)  
363 //       << " ; hasMother = " << hasMother << endl;
364
365     Bool_t hasDaughter = (iparticle->GetNDaughters() > 0);
366
367    // cout << "AliGenUHKM::Generate() n.daughters = " << iparticle->GetNDaughters() 
368     //<< " ; hasDaughter = " << hasDaughter << endl;
369
370
371     if(hasDaughter) {
372       //      cout << "AliGenUHKM::Generate() decayed particle (not trackable)" << endl;
373       // This particle has decayed
374       // It will not be tracked
375       // Add it only once with coordiinates not
376       // smeared with primary vertex position
377       Float_t p[3] = {p[0] = iparticle->Px(),
378                       p[1] = iparticle->Py(),
379                       p[2] = iparticle->Pz()};
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       v[0] = iparticle->Vx();
383       v[1] = iparticle->Vy();
384       v[2] = iparticle->Vz();
385       Float_t time = iparticle->T();
386
387       Int_t type    = iparticle->GetStatusCode(); //j1/h0
388
389       Int_t imo = -1;
390       if(hasMother) {
391         imo = iparticle->GetFirstMother(); //index of mother particle in fParticles
392       } // if has mother
393       Bool_t trackFlag = (!hasDaughter);   // tFlag is kFALSE --> do not track the particle
394
395 //      printf("Pushing Track %d with status %d mother %d\n", kf, trackFlag, imo>=0?idsOnStack[imo]:imo);
396       PushTrack(trackFlag, (imo>=0 ? idsOnStack[imo+1] : imo), kf,
397                 p[0], p[1], p[2], energy,
398                 v[0], v[1], v[2], time,
399                 polar[0], polar[1], polar[2],
400                 (hasMother ? kPDecay : kPNoProcess), nt);
401     //  cout << "AliGenUHKM::Generate() pushed on stack with stack index = " << nt 
402 //         << "; mother index on stack = " << (imo>=0 ? idsOnStack[imo+1] : imo) << endl;
403       idsOnStack[i] = nt;
404       fNprimaries++;
405       KeepTrack(nt);
406     }
407     else {
408       //      cout << "AliGenUHKM::Generate() final particle --> push it twice on the stack" << endl;
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 ccordiantes and vertex smearing
415       //   this one will be tracked
416       Float_t p[3] = {p[0] = iparticle->Px(),
417                       p[1] = iparticle->Py(),
418                       p[2] = iparticle->Pz()};
419       mass = TDatabasePDG::Instance()->GetParticle(kf)->Mass();
420       energy = sqrt(mass*mass + p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
421       v[0] = iparticle->Vx();
422       v[1] = iparticle->Vy();
423       v[2] = iparticle->Vz();
424
425       Int_t type    = iparticle->GetStatusCode(); //j1/h0
426       Int_t coeffT=1;
427       if(type==1)coeffT=-1; //to separate particles from jets
428
429       Int_t imo = -1;
430       
431       if(hasMother) {
432         imo = iparticle->GetFirstMother();
433       } // if has mother
434       Bool_t trackFlag = (hasDaughter);  // tFlag = kFALSE --> do not track this one, its for femtoscopy
435        
436       PushTrack(trackFlag, (imo>=0 ? idsOnStack[imo+1] : imo), kf,
437                 p[0], p[1], p[2], energy,
438                 v[0], v[1], v[2], (iparticle->T())*coeffT,
439                 polar[0], polar[1], polar[2],
440                 hasMother ? kPDecay:kPNoProcess, nt);
441      // cout << "AliGenUHKM::Generate() pushed on stack with stack index = " << nt
442     //       << "; mother index on stack = " << (imo>=0 ? idsOnStack[imo+1] : imo) << endl;
443
444       idsOnStack[i] = nt;
445       fNprimaries++;
446       KeepTrack(nt);
447
448       origin[0] = origin0[0]+v[0];
449       origin[1] = origin0[1]+v[1];
450       origin[2] = origin0[2]+v[2];
451       imo = nt;
452       
453       trackFlag = (!hasDaughter);    // tFlag = kTRUE --> track this one
454       //cout << "AliGenUHKM::Generate() trackFlag = " << trackFlag << endl;
455
456       PushTrack(trackFlag, imo, kf,
457                 p[0], p[1], p[2], energy,
458                 origin[0], origin[1], origin[2], iparticle->T(),
459                 polar[0], polar[1], polar[2],
460                 hasMother ? kPDecay:kPNoProcess, nt);
461      // cout << "AliGenUHKM::Generate() pushed on stack with stack index = " << nt
462     //       << "; mother index on stack = " << imo << endl;
463       fNprimaries++;
464       KeepTrack(nt);
465     }
466   }
467
468   SetHighWaterMark(fNprimaries);
469
470   TArrayF eventVertex;
471   eventVertex.Set(3);
472   eventVertex[0] = origin0[0];
473   eventVertex[1] = origin0[1];
474   eventVertex[2] = origin0[2];
475
476   // Builds the event header, to be called after each event
477   AliGenEventHeader* header = new AliGenHijingEventHeader("UHKM");
478
479   ((AliGenHijingEventHeader*) header)->SetNProduced(fNprimaries);
480   ((AliGenHijingEventHeader*) header)->SetPrimaryVertex(eventVertex);
481   ((AliGenHijingEventHeader*) header)->SetImpactParameter(0.0);
482   ((AliGenHijingEventHeader*) header)->SetTotalEnergy(0.0);
483   ((AliGenHijingEventHeader*) header)->SetHardScatters(0);
484   ((AliGenHijingEventHeader*) header)->SetParticipants(0, 0);
485   ((AliGenHijingEventHeader*) header)->SetCollisions(0, 0, 0, 0);
486   ((AliGenHijingEventHeader*) header)->SetSpectators(0, 0, 0, 0);
487   ((AliGenHijingEventHeader*) header)->SetReactionPlaneAngle(0);//evrot);
488
489   header->SetPrimaryVertex(fVertex);
490   AddHeader(header);
491   fCollisionGeometry = (AliGenHijingEventHeader*)  header;
492
493   delete idsOnStack;
494
495   //  gAlice->SetGenEventHeader(header);
496
497   printf(" Finish Generate .. %d ..\n",nt);
498   cout << "AliGenUHKM::Generate() OUT" << endl;
499 }
500
501 void AliGenUHKM::Copy(TObject &) const
502 {
503   Fatal("Copy","Not implemented!\n");
504 }
505
506 void AliGenUHKM::SetAllParameters() {
507   cout << "AliGenUHKM::SetAllParameters() IN" << endl;
508
509   fUHKMgen->SetEcms(fHydjetParams.fSqrtS);
510   fUHKMgen->SetBmin(fHydjetParams.fBmin);
511   fUHKMgen->SetBmax(fHydjetParams.fBmax);
512   fUHKMgen->SetAw(fHydjetParams.fAw);
513   fUHKMgen->SetSeed(fHydjetParams.fSeed);
514
515   fUHKMgen->SetChFrzTemperature(fHydjetParams.fT);
516   fUHKMgen->SetMuB(fHydjetParams.fMuB);
517   fUHKMgen->SetMuS(fHydjetParams.fMuS);
518   fUHKMgen->SetMuQ(fHydjetParams.fMuI3);
519   fUHKMgen->SetTauB(fHydjetParams.fTau);
520   fUHKMgen->SetThFrzTemperature(fHydjetParams.fThFO);
521   fUHKMgen->SetMuPionThermal(fHydjetParams.fMu_th_pip);
522
523   fUHKMgen->SetSigmaTau(fHydjetParams.fSigmaTau);
524   fUHKMgen->SetRmaxB(fHydjetParams.fR);
525   fUHKMgen->SetYlMax(fHydjetParams.fYlmax);
526   fUHKMgen->SetEtaRMax(fHydjetParams.fUmax);
527   fUHKMgen->SetMomAsymmPar(fHydjetParams.fDelta);
528   fUHKMgen->SetCoordAsymmPar(fHydjetParams.fEpsilon);
529
530   fUHKMgen->SetGammaS(fHydjetParams.fCorrS);
531   // fUHKMgen->SetHdec(fHydjetParams.fTime);
532   fUHKMgen->SetEtaType(fHydjetParams.fEtaType);
533   fUHKMgen->SetFlagWeakDecay(fHydjetParams.fWeakDecay);
534
535   //PYQUEN parameters
536
537   fUHKMgen->SetPyquenNhsel(fHydjetParams.fNhsel);
538   fUHKMgen->SetPyquenShad(fHydjetParams.fIshad);
539   fUHKMgen->SetPyquenPtmin(fHydjetParams.fPtmin);
540   fUHKMgen->SetPyquenT0(fHydjetParams.fT0);
541   fUHKMgen->SetPyquenTau0(fHydjetParams.fTau0);
542   fUHKMgen->SetPyquenNf(fHydjetParams.fNf);
543   fUHKMgen->SetPyquenIenglu(fHydjetParams.fIenglu);
544   fUHKMgen->SetPyquenIanglu(fHydjetParams.fIanglu);
545
546   fUHKMgen->SetPDGParticleFile(fParticleFilename);
547   fUHKMgen->SetPDGDecayFile(fDecayFilename);
548   for(Int_t i=0; i<fStableFlagged; i++) fUHKMgen->SetPDGParticleStable(fStableFlagPDG[i], fStableFlagStatus[i]);
549   fUHKMgen->SetUseCharmParticles(fUseCharmParticles);
550   fUHKMgen->SetMinimumWidth(fMinWidth);
551   fUHKMgen->SetMaximumWidth(fMaxWidth);
552   fUHKMgen->SetMinimumMass(fMinMass);
553   fUHKMgen->SetMaximumMass(fMaxMass);
554
555  cout<<" Print all parameters "<<endl;
556  cout<<" SqrtS = "<<fHydjetParams.fSqrtS<<endl;
557  cout<<" Bmin  = "<< fHydjetParams.fBmin<<endl;
558  cout<<" Bmax= "<<fHydjetParams.fBmax<<endl;
559  cout<<" Aw= "<<fHydjetParams.fAw<<endl;
560  cout<<" Seed= "<<fHydjetParams.fSeed<<endl;
561
562  cout<<" ---Stat-model parameters----------- "<<endl;
563
564  cout<<" ChFrzTemperature= "<<fHydjetParams.fT<<endl;
565  cout<<" MuB= "<<fHydjetParams.fMuB<<endl;
566  cout<<" MuS= "<<fHydjetParams.fMuS<<endl;
567  cout<<" MuQ= "<<fHydjetParams.fMuI3<<endl;
568  cout<<" TauB= "<<fHydjetParams.fTau<<endl;
569  cout<<" ThFrzTemperature= "<<fHydjetParams.fThFO<<endl;
570  cout<<" MuPionThermal= "<<fHydjetParams.fMu_th_pip<<endl;
571
572 cout<<"-----Volume parameters -------------- "<<endl;
573
574  cout<<" SigmaTau= "<<fHydjetParams.fSigmaTau<<endl;
575  cout<<" RmaxB= "<<fHydjetParams.fR<<endl;
576  cout<<" YlMax= "<<fHydjetParams.fYlmax<<endl;
577  cout<<" EtaRMax= "<<fHydjetParams.fUmax<<endl;
578  cout<<" MomAsymmPar= "<<fHydjetParams.fDelta<<endl;
579  cout<<" CoordAsymmPar= "<<fHydjetParams.fEpsilon<<endl;
580
581 cout<<" --------Flags------ "<<endl;
582
583  cout<<" GammaS= "<<fHydjetParams.fCorrS<<endl;
584   // fUHKMgen->SetHdec(fHydjetParams.fTime<<endl;
585  cout<<" EtaType= "<<fHydjetParams.fEtaType<<endl;
586  cout<<" FlagWeakDecay= "<<fHydjetParams.fWeakDecay<<endl;
587
588   cout<<"----PYQUEN parameters---"<<endl;
589
590   cout<<" Nhsel= "<<fHydjetParams.fNhsel<<endl;
591   cout<<" Shad= "<<fHydjetParams.fIshad<<endl;
592   cout<<" Ptmin= "<<fHydjetParams.fPtmin<<endl;
593   cout<<" T0= "<<fHydjetParams.fT0<<endl;
594   cout<<" Tau0= "<<fHydjetParams.fTau0<<endl;
595   cout<<" Nf= "<<fHydjetParams.fNf<<endl;
596   cout<<" Ienglu= "<<fHydjetParams.fIenglu<<endl;
597   cout<<" Ianglu= "<<fHydjetParams.fIanglu<<endl;
598
599   cout<<"----PDG table parameters---"<<endl;
600   
601   cout<<" UseCharmParticles= "<<fUseCharmParticles<<endl;
602   cout<<" MinimumWidth= "<<fMinWidth<<endl;
603   cout<<" MaximumWidth= "<<fMaxWidth<<endl;
604   cout<<" MinimumMass= "<<fMinMass<<endl;
605   cout<<" MaximumMass= "<<fMaxMass<<endl;
606
607
608
609   cout << "AliGenUHKM::SetAllParameters() OUT" << endl;
610 }
611
612 // add the additional PDG codes from UHKM(SHARE table) to ROOT's table
613 void AliGenUHKM::CheckPDGTable() {
614   cout << "AliGenUHKM::CheckPDGTable()   IN" << endl;
615   //TDabasePDG *rootPDG  = TDatabasePDG::Instance();         // ROOT's PDG table
616   DatabasePDG *uhkmPDG = fUHKMgen->PDGInfo();              // UHKM's PDG table
617   TParticlePDG *rootTestParticle;
618   ParticlePDG *uhkmTestParticle;
619
620   cout << "particles with good status in UHKM table = " << uhkmPDG->GetNParticles() << endl;
621   // loop over all particles in the SHARE table
622   for(Int_t i=0; i<uhkmPDG->GetNParticles(); i++) {
623     cout << "particle #" << i << " ================" << endl;
624     // get a particle specie
625     uhkmTestParticle = uhkmPDG->GetPDGParticleByIndex(i);
626     cout << "PDG = " << uhkmTestParticle->GetPDG() << endl;
627     // check if this code exists in ROOT's table
628     rootTestParticle = TDatabasePDG::Instance()->GetParticle(uhkmTestParticle->GetPDG());
629     if(!rootTestParticle) {    // if not then add it to the ROOT's PDG database
630       
631       TDatabasePDG::Instance()->AddParticle(uhkmTestParticle->GetName(), uhkmTestParticle->GetName(), 
632                                             uhkmTestParticle->GetMass(), uhkmTestParticle->GetElectricCharge(),
633                                             (uhkmTestParticle->GetWidth()<1e-10 ? kTRUE : kFALSE),
634                                             uhkmTestParticle->GetWidth(), 
635                                             (Int_t(uhkmTestParticle->GetBaryonNumber())==0 ? "meson" : "baryon"),
636                                             uhkmTestParticle->GetPDG());    
637      cout << "Not found in ROOT's PDG database --> added now" << endl;
638     if(uhkmTestParticle->GetWidth()<1e-10) cout<<uhkmTestParticle->GetPDG()<<"its mass "<< 
639     TDatabasePDG::Instance()->GetParticle(uhkmTestParticle->GetPDG())->Mass()<<
640     TDatabasePDG::Instance()->GetParticle(uhkmTestParticle->GetPDG())->Width()<<endl;  
641     }
642     else
643       cout << "Found in ROOT's PDG database --> not added" << endl;
644   }  // end for
645   cout << "AliGenUHKM::CheckPDGTable()   OUT" << endl;
646 }