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