TUHKMgen/TUHKMgen.cxx

   1 /////////////////////////////////////////////////////////////////////////////
   2 // TUHKM is an interface to
   3 // UHKM 3.0                                                                //
   4 // ( only the HYDJET++ part is implemented)                                //
   5 // temporary link:                                                         //
   6 // http://lav01.sinp.msu.ru/~igor/hydjet++/hydjet++.txt                    //
   7 // The main UHKM options are accessable  through this interface.           //
   8 // Uses the TUHKMgen implementation of TGenerator.                         //
   9 // Author of the first implementation: Sergey Zaporozhets                  //
  10 // (zaporozh@sunhe.jinr.ru)                                                //
  11 // Futhers modifications were made by                                      //
  12 // Ionut Cristian Arsene (i.c.arsene@fys.uio.no)                           //
  13 // & Malinina Liudmila(malinina@lav01.sinp.msu.ru) using as an example     //
  14 //  AliGenTherminator.cxx created by Adam Kisiel                           //
  15 //                                                                         //
  16 ////////////////////////////////////////////////////////////////////////////
  17
  18 #ifndef TUHKMGEN_H
  19 #include "TUHKMgen.h"
  20 #endif
  21 #include "TObjArray.h"
  22 #include "TParticle.h"
  23 #include "Particle.h"
  24 #include "TROOT.h"
  25 #include "TError.h"
  26 #include "TClonesArray.h"
  27 #include "TSystem.h"
  28
  29 #include <string>
  30 #include <iostream>
  31 using namespace std;
  32
  33 //class TObjArray;
  34
  35 ClassImp(TUHKMgen)
  36
  37 TUHKMgen::TUHKMgen() :
  38   TGenerator("UHKM","UHKM"),
  39   fInitialState(0x0),
  40   fAllocator(),
  41   fSourceList(),
  42   fSecondariesList(),
  43   fNPprim(0),
  44   fNPsec(0),
  45   fHydjetParams(),
  46   fStableFlagged(0)
  47   //  fUseCharmParticles(kFALSE),
  48   //  fMinWidth(0.0),
  49   //  fMaxWidth(10.0),
  50   //  fMinMass(0.0001),
  51   //  fMaxMass(10.0)
  52 {
  53   // default constructor setting reasonable defaults for initial parameters (central Pb+Pb at 5.5 TeV)
  54
  55   //  cout << "TUHKMgen::TUHKMgen() IN" << endl;
  56
  57   // Set reasonable default values for LHC
  58
  59   fHydjetParams.fSqrtS=5500; //LHC
  60   fHydjetParams.fAw=207;//Au-Au
  61   fHydjetParams.fIfb    = 1;      // flag of type of centrality generation (=0 is fixed by fBfix, not 0
  62                                    //impact parameter is generated in each event between fBfmin
  63                                    //and fBmax according with Glauber model (f-la 30)
  64   fHydjetParams.fBfix=0.;
  65   fHydjetParams.fBmin=0.;
  66   fHydjetParams.fBmax=0.5; //0-5% centrality
  67   fHydjetParams.fT = 0.170;
  68   fHydjetParams.fMuB = 0.0;
  69   fHydjetParams.fMuS = 0.0;
  70   fHydjetParams.fMuI3 = 0.0;
  71   fHydjetParams.fThFO = 0.130;
  72   fHydjetParams.fMu_th_pip = 0.0;
  73   fHydjetParams.fSeed=0;
  74   fHydjetParams.fTau=10.;
  75   fHydjetParams.fSigmaTau=3.;
  76   fHydjetParams.fR=11.;
  77   fHydjetParams.fYlmax=4.0;
  78   fHydjetParams.fUmax=1.1;
  79   fHydjetParams.fDelta=0.;
  80   fHydjetParams.fEpsilon=0.;
  81   fHydjetParams.fWeakDecay=0; //>=0 on ,-1 off
  82   fHydjetParams.fEtaType=1;//gaus
  83   fHydjetParams.fCorrS=1.;
  84   fHydjetParams.fNhsel=2;
  85   fHydjetParams.fIshad=1;
  86   fHydjetParams.fPtmin=7.0;
  87   fHydjetParams.fT0=0.8;
  88   fHydjetParams.fTau0=0.1;
  89   fHydjetParams.fNf=0;
  90   fHydjetParams.fIenglu=0;
  91   fHydjetParams.fIanglu=0;
  92
  93
  94   // Set reasonable default values for RHIC
  95 /*
  96   fHydjetParams.fSqrtS=200; //RHIC
  97   fHydjetParams.fAw=197;//Au-Au
  98   fHydjetParams.fIfb    = 1;
  99   fHydjetParams.fBfix=0.;
 100   fHydjetParams.fBmin=0.;
 101   fHydjetParams.fBmax=0.5; //0-5% centrality
 102   fHydjetParams.fT = 0.165;
 103   fHydjetParams.fMuB = 0.0285;
 104   fHydjetParams.fMuS = 0.007;
 105   fHydjetParams.fMuI3 = -0.001;
 106   fHydjetParams.fThFO = 0.100;
 107   fHydjetParams.fMu_th_pip = 0.053;
 108   fHydjetParams.fSeed=0;
 109   fHydjetParams.fTau=8.;
 110   fHydjetParams.fSigmaTau=2.;
 111   fHydjetParams.fR=10.;
 112   fHydjetParams.fYlmax=3.3;
 113   fHydjetParams.fUmax=1.1;
 114   fHydjetParams.fDelta=0.;
 115   fHydjetParams.fEpsilon=0.;
 116   fHydjetParams.fWeakDecay=0;//>=0 on ,-1 off
 117   fHydjetParams.fEtaType=1;//gaus
 118   fHydjetParams.fCorrS=1.;
 119   fHydjetParams.fNhsel=2;
 120   fHydjetParams.fIshad=1;
 121   fHydjetParams.fPtmin=3.4;
 122   fHydjetParams.fT0=0.3;
 123   fHydjetParams.fTau0=0.4;
 124   fHydjetParams.fNf=2;
 125   fHydjetParams.fIenglu=0;
 126   fHydjetParams.fIanglu=0;
 127 */
 128
 129   strcpy(fParticleFilename, Form("%s/TUHKMgen/UHKM/particles.data", gSystem->Getenv("ALICE_ROOT")));
 130   strcpy(fDecayFilename, Form("%s/TUHKMgen/UHKM/tabledecay.txt", gSystem->Getenv("ALICE_ROOT")));
 131   for(Int_t i=0; i<500; i++) {
 132     fStableFlagPDG[i] = 0;
 133     fStableFlagStatus[i] = kFALSE;
 134   }
 135   fStableFlagged = 0;
 136   //  cout << "TUHKMgen::TUHKMgen() OUT" << endl;
 137 }
 138
 139 //______________________________________________________________________________
 140 TUHKMgen::~TUHKMgen()
 141 {
 142   // destructor, deletes the InitialStateHydjet object
 143   //  cout << "TUHKMgen::~TUHKMgen() IN" << endl;
 144   if(fInitialState)
 145     delete fInitialState;
 146   //  cout << "TUHKMgen::~TUHKMgen() OUT" << endl;
 147 }
 148
 149 void TUHKMgen::SetAllParametersRHIC()
 150 {
 151   // Set reasonable input parameters for 0-5% central Au+Au collisions
 152   // at 200 GeV at RHIC
 153   SetEcms(200.0);                // RHIC top energy
 154   SetAw(197);                    // Au+Au
 155   SetIfb(1);
 156   SetBfix(0.);
 157   SetBmin(0.0);                  // 0%
 158   SetBmax(0.5);                  // 5%
 159   SetChFrzTemperature(0.165);    // T_ch = 165 MeV
 160   SetMuB(0.0285);                // mu_B = 28.5 MeV
 161   SetMuS(0.007);                 // mu_S = 7 MeV
 162   SetMuQ(-0.001);                // mu_Q = -1 MeV
 163   SetThFrzTemperature(0.100);    // T_th = 100 MeV
 164   SetMuPionThermal(0.053);       // mu_th_pion = 53 MeV
 165   SetSeed(0);                    // use UNIX time
 166   SetTauB(8.0);                  // tau = 8 fm/c
 167   SetSigmaTau(2.0);              // sigma_tau = 2 fm/c
 168   SetRmaxB(10.0);                // fR = 10 fm
 169   SetYlMax(3.3);                 // fYmax = 3.3
 170   SetEtaRMax(1.1);               // Umax = 1.1
 171   SetMomAsymmPar(0.0);           // delta = 0.0
 172   SetCoordAsymmPar(0.0);         // epsilon = 0.0
 173   SetFlagWeakDecay(0.0);           // weak decay on (<0 off !!!)
 174   SetEtaType(1);                 // gaus distributed with fYmax dispersion (0 means boost invariant)
 175   SetGammaS(1.0);                // gammaS = 1.0 (no strangeness canonical suppresion)
 176   SetPyquenNhsel(2);             // hydro on, jets on, jet quenching on
 177   SetPyquenShad(1);              // shadowing on (0 off)
 178   SetPyquenPtmin(3.4);           // ptmin = 3.4 GeV/c
 179   SetPyquenT0(0.3);              // T0 = 300 MeV
 180   SetPyquenTau0(0.4);            // tau0 = 0.4 fm/c
 181   SetPyquenNf(2);                // 2 flavours
 182   SetPyquenIenglu(0);            // radiative and collisional energy loss
 183   SetPyquenIanglu(0);            // small gluon angular distribution
 184 }
 185
 186 void TUHKMgen::SetAllParametersLHC()
 187 {
 188   // Set reasonable input parameters for 0-5% Pb+Pb collisions at 5.5 TeV at LHC
 189   SetEcms(5500.0);               // LHC
 190   SetAw(207);                    // Pb+Pb
 191   SetIfb(1);
 192   SetBfix(0.);                  // 0
 193   SetBmin(0.0);                  // 0%
 194   SetBmax(0.5);                  // 5%
 195   SetChFrzTemperature(0.170);    // T_ch = 170 MeV
 196   SetMuB(0.0);                   // mu_B = 0 MeV
 197   SetMuS(0.0);                   // mu_S = 0 MeV
 198   SetMuQ(0.0);                   // mu_Q = 0 MeV
 199   SetThFrzTemperature(0.130);    // T_th = 130 MeV
 200   SetMuPionThermal(0.0);         // mu_th_pion = 0 MeV
 201   SetSeed(0);                    // use UNIX time
 202   SetTauB(10.0);                 // tau = 10 fm/c
 203   SetSigmaTau(3.0);              // sigma_tau = 3 fm/c
 204   SetRmaxB(11.0);                // fR = 11 fm
 205   SetYlMax(4.0);                 // fYmax = 4.0
 206   SetEtaRMax(1.1);               // Umax = 1.1
 207   SetMomAsymmPar(0.0);           // delta = 0.0
 208   SetCoordAsymmPar(0.0);         // epsilon = 0.0
 209   SetFlagWeakDecay(0);           // weak decay on (<0 off !!!)
 210   SetEtaType(1);                 // gaus distributed with fYmax dispersion (0 means boost invariant)
 211   SetGammaS(1.0);                // gammaS = 1.0 (no strangeness canonical suppresion)
 212   SetPyquenNhsel(2);             // hydro on, jets on, jet quenching on
 213   SetPyquenShad(1);              // shadowing on (0 off)
 214   SetPyquenPtmin(7.0);           // ptmin = 7.0 GeV/c
 215   SetPyquenT0(0.8);              // T0 = 800 MeV
 216   SetPyquenTau0(0.1);            // tau0 = 0.4 fm/c
 217   SetPyquenNf(0);                // 0 flavours
 218   SetPyquenIenglu(0);            // radiative and collisional energy loss
 219   SetPyquenIanglu(0);            // small gluon angular distribution
 220 }
 221
 222 TObjArray* TUHKMgen::ImportParticles(const Option_t *)
 223 {
 224   // Function overloading the TGenerator::ImportParticles() member function.
 225   // The particles from the local particle list (fSecondariesList) are
 226   // forwarded to the TGenerator::fParticles
 227   //  cout << "TObjArray* TUHKMgen::ImportParticles(Option_t) IN" << endl;
 228   fParticles->Clear();
 229   //  cout << "TObjArray* TUHKMgen::ImportParticles(Option_t): UHKM stack contains " << fNPsec << " particles." << endl;
 230   Int_t nump = 0;
 231   LPIT_t it,e;
 232
 233   for(it = fSecondariesList.begin(), e = fSecondariesList.end(); it != e; ++it) {
 234     TVector3 pos(it->Pos().Vect());
 235     TVector3 mom(it->Mom().Vect());
 236     Float_t m1 = it->TableMass();
 237     Int_t im1 = it->GetMother();
 238     Int_t im2 = -1;
 239     Int_t id1 = -1;
 240     Int_t id2 = -1;
 241
 242     //    Int_t nd = it->GetNDaughters();
 243
 244     Int_t type = it->GetType();  // 0-hydro, 1-jets
 245
 246     if (im1> -1) {
 247       TParticle *mother = (TParticle*) (fParticles->UncheckedAt(im1));
 248       mother->SetLastDaughter(nump);
 249       if (mother->GetFirstDaughter()==-1)
 250         mother->SetFirstDaughter(nump);
 251     }
 252
 253     nump++;
 254
 255     TParticle* p = new TParticle(it->Encoding(), type,                          //pdg,stat
 256                                  im1, im2, id1, id2,                                            //m1,m2,d1,d2
 257                                  mom[0], mom[1], mom[2], TMath::Sqrt(mom.Mag2() + m1 * m1),     //px,py,pz,e
 258                                  pos[0]*1.e-13, pos[1]*1.e-13, pos[2]*1.e-13,
 259                                  it->T()*1.e-13/3e10);                          //x,y,z,t
 260
 261      p->SetUniqueID(nump);
 262      fParticles->Add(p);
 263   } //end for
 264
 265   fAllocator.FreeList(fSourceList);
 266   fAllocator.FreeList(fSecondariesList);
 267   //  cout << "TObjArray* TUHKMgen::ImportParticles(Option_t) OUT" << endl;
 268   return fParticles;
 269 }
 270
 271 Int_t TUHKMgen::ImportParticles(TClonesArray *particles, const Option_t* option)
 272 {
 273   // Function overloading the TGenerator::ImportParticles() member function.
 274   // The particles from the local particle list (fSecondariesList) are
 275   // forwarded to the TGenerator::fParticles
 276   //  cout << "Int_t TUHKMgen::ImportParticles(TClonesArray*, Option_t) IN" << endl;
 277   if(particles==0) return 0;
 278   TClonesArray &particlesR=*particles;
 279   particlesR.Clear();
 280
 281
 282   Int_t numprim,numsec;  numprim=numsec=0;
 283   Int_t nump = 0;
 284   LPIT_t it,e;
 285
 286   //  cout << "TUHKMgen::ImportParticles() option(All or Sec) = " << option << endl;
 287   for(it = fSecondariesList.begin(), e = fSecondariesList.end(); it != e; ++it) {
 288     //!!! for(Int_t pp=0;pp<100;pp++) {
 289         //      cout << "TUHKMgen::ImportParticles() import particle pdg(" << it->Encoding() << ")" << endl;
 290     TVector3 pos(it->Pos().Vect());  //
 291     //!!! TVector3 pos(0.0, 0.1, 0.2);
 292     TVector3 mom(it->Mom().Vect());  //
 293     //!!! TVector3 mom(-0.1, 0.1, 0.2);
 294     Float_t m1 = it->TableMass();  //
 295     //!!!    Float_t m1 = 0.139;
 296     Int_t im1 = it->GetMother();  //
 297     //!!!    Int_t im1 = (pp<50? -1 : pp-50);
 298     //    Int_t nd = 0;
 299     Int_t im2 = -1;
 300     Int_t id1 = -1;
 301     Int_t id2 = -1;
 302
 303     //    nd = it->GetNDaughters();
 304
 305     Int_t type = it->GetType();  //  // 0-hydro, 1-jets
 306     //!!!    Int_t type = 0;
 307
 308     if (im1> -1) {
 309       TParticle *mother = (TParticle*) (particlesR.UncheckedAt(im1+1));
 310       mother->SetLastDaughter(nump);
 311       if (mother->GetFirstDaughter()==-1)
 312         mother->SetFirstDaughter(nump);
 313     }
 314
 315     nump++;
 316
 317     new (particlesR[nump]) TParticle(it->Encoding(), type,                                              //pdg,stat
 318                                      im1, im2, id1, id2,                                                //m1,m2,d1,d2
 319                                      mom[0], mom[1], mom[2], TMath::Sqrt(mom.Mag2() + m1 * m1), //px,py,pz,e
 320                                      pos[0]*1.e-13, pos[1]*1.e-13, pos[2]*1.e-13,
 321                                      it->T()*1.e-13/3e10);                      //x,y,z,t
 322
 323     /*!!! new (particlesR[nump]) TParticle(211, type,                                         //pdg,stat
 324                                      im1, im2, id1, id2,                                                //m1,m2,d1,d2
 325                                      mom[0], mom[1], mom[2], TMath::Sqrt(mom.Mag2() + m1 * m1), //px,py,pz,e
 326                                      pos[0]*1.e-13, pos[1]*1.e-13, pos[2]*1.e-13,
 327                                      0.0*1.e-13/3e10);                       //x,y,z,t
 328     !!!*/
 329
 330     particlesR[nump]->SetUniqueID(nump);
 331     numsec++;
 332   }//end for
 333
 334   fAllocator.FreeList(fSourceList);  //
 335   fAllocator.FreeList(fSecondariesList);  //
 336   //  printf("Scan and add prim %d sec %d and all %d particles\n",
 337   //     numprim,numsec,nump);
 338   //  cout << "Int_t TUHKMgen::ImportParticles(TClonesArray*, Option_t) OUT" << endl;
 339   return nump;
 340 }
 341
 342 //______________________________________________________________________________
 343 void TUHKMgen::Initialize()
 344 {
 345   // Function overloading the TGenerator::Initialize() member function.
 346   // The Monte-Carlo model is initialized (input parameters are transmitted,
 347   // particle list and decay channels are loaded, average multiplicities are calculated, etc.)
 348   //  cout << "TUHKMgen::Initialize() IN" << endl;
 349   fInitialState = new InitialStateHydjet();  //
 350   SetAllParameters();  //
 351   fInitialState->LoadPDGInfo();  //
 352   // set the stable flags
 353   for(Int_t i=0; i<fStableFlagged; i++) //
 354     fInitialState->SetPDGParticleStable(fStableFlagPDG[i], fStableFlagStatus[i]); //
 355
 356   if(!fInitialState->MultIni()) //
 357     Error("TUHKMgen::Initialize()", "Bad status return from MultIni(). Check it out!! \n");  //
 358   //  cout << "TUHKMgen::Initialize() OUT" << endl;
 359 }
 360
 361 void TUHKMgen::Print(const Option_t*) const
 362 {
 363   cout << "TUHKMgen::Print() method not implemented yet!!" << endl;
 364 }
 365
 366 void TUHKMgen::GenerateEvent()
 367 {
 368   // Member function overloading the TGenerator::GenerateEvent()
 369   // The HYDJET++ model is run and the particle lists (fSourceList and fSecondariesList) are filled
 370   //  cout << "TUHKMgen::GenerateEvent() IN" << endl;
 371  // cout << "TUHKMgen::GenerateEvent() IN fSourceList " <<endl;
 372   //  cout<< " fInitialStrate "<<fInitialState<<endl;  //
 373   //" fSourceList "<< fSourceList<<" fAllocator "<< fAllocator <<endl;
 374
 375   // Generate the initial particles
 376   fInitialState->Initialize(fSourceList, fAllocator);  //
 377   //  cout << "TUHKMgen::fInitialState->Initialize" << endl;
 378
 379   if(fSourceList.empty()) //
 380     Error("TUHKMgen::GenerateEvent()", "Source particle list empty after fireball initialization!! \n");  //
 381
 382   // Run the decays
 383
 384   //  cout << "after Initilize: get_time, weak_decay_limit" <<fInitialState->GetTime()<<fInitialState->GetWeakDecayLimit()<< endl;   //
 385   if(fInitialState->GetTime() >= 0.)  //
 386     fInitialState->Evolve(fSourceList, fSecondariesList, fAllocator, fInitialState->GetWeakDecayLimit());  //
 387   //  cout << "TUHKMgen::GenerateEvent() OUT" << endl;
 388 }
 389
 390 void TUHKMgen::SetAllParameters() {
 391   // forward all model input parameters to the InitialStateHydjet object
 392   // which will handle all the Monte-Carlo simulation using HYDJET++ model
 393
 394   //  cout << "TUHKMgen::SetAllParameters() IN" << endl;
 395
 396   fInitialState->fParams.fSqrtS = fHydjetParams.fSqrtS;
 397   fInitialState->fParams.fAw = fHydjetParams.fAw;
 398   fInitialState->fParams.fIfb = fHydjetParams.fIfb ;
 399   fInitialState->fParams.fBfix = fHydjetParams.fBfix;
 400   fInitialState->fParams.fBmin = fHydjetParams.fBmin;
 401   fInitialState->fParams.fBmax = fHydjetParams.fBmax;
 402   fInitialState->fParams.fSeed = fHydjetParams.fSeed;
 403   fInitialState->fParams.fT = fHydjetParams.fT;
 404   fInitialState->fParams.fMuB = fHydjetParams.fMuB;
 405   fInitialState->fParams.fMuS = fHydjetParams.fMuS;
 406   fInitialState->fParams.fMuI3 = fHydjetParams.fMuI3;
 407   fInitialState->fParams.fThFO = fHydjetParams.fThFO;
 408   fInitialState->fParams.fMu_th_pip = fHydjetParams.fMu_th_pip;
 409
 410   fInitialState->fParams.fTau = fHydjetParams.fTau;
 411   fInitialState->fParams.fSigmaTau = fHydjetParams.fSigmaTau;
 412   fInitialState->fParams.fR = fHydjetParams.fR;
 413   fInitialState->fParams.fYlmax = fHydjetParams.fYlmax;
 414   fInitialState->fParams.fUmax = fHydjetParams.fUmax;
 415   fInitialState->fParams.fDelta = fHydjetParams.fDelta;
 416   fInitialState->fParams.fEpsilon = fHydjetParams.fEpsilon;
 417
 418  // fInitialState->fParams.fTime = fHydjetParams.fTime;
 419   fInitialState->fParams.fWeakDecay = fHydjetParams.fWeakDecay;
 420   fInitialState->fParams.fCorrS = fHydjetParams.fCorrS;
 421   fInitialState->fParams.fEtaType = fHydjetParams.fEtaType;
 422
 423   fInitialState->fParams.fPtmin = fHydjetParams.fPtmin;
 424   fInitialState->fParams.fNhsel = fHydjetParams.fNhsel;
 425   fInitialState->fParams.fIshad = fHydjetParams.fIshad;
 426   fInitialState->fParams.fT0 = fHydjetParams.fT0;
 427   fInitialState->fParams.fTau0 = fHydjetParams.fTau0;
 428   fInitialState->fParams.fNf = fHydjetParams.fNf;
 429   fInitialState->fParams.fIenglu = fHydjetParams.fIenglu;
 430   fInitialState->fParams.fIanglu = fHydjetParams.fIanglu;
 431
 432   fInitialState->SetPDGParticleFilename(fParticleFilename);
 433   fInitialState->SetPDGDecayFilename(fDecayFilename);
 434   //  fInitialState->SetUseCharmParticles(fUseCharmParticles);
 435   //  fInitialState->SetWidthRange(fMinWidth, fMaxWidth);
 436   //  fInitialState->SetMassRange(fMinMass, fMaxMass);
 437   //  for(Int_t i=0; i<fStableFlagged; i++)
 438   //    fInitialState->SetPDGParticleStable(fStableFlagPDG[i], fStableFlagStatus[i]);
 439   //  cout << "TUHKMgen::SetAllParameters() OUT" << endl;
 440 }
 441