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 }