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