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   fNPprim(0),
  41   fNPsec(0),
  42   fHydjetParams(),
  43   fStableFlagged(0)
  44 {
  45   // default constructor setting reasonable defaults for initial parameters (central Pb+Pb at 5.5 TeV)
  46
  47   ParticleAllocator fAllocator;
  48   List_t fSecondariesList;
  49
  50   // Set reasonable default values for LHC
  51
  52   fHydjetParams.fSqrtS=5500; //LHC
  53   fHydjetParams.fAw=207;//Au-Au
  54   fHydjetParams.fIfb    = 1;      // flag of type of centrality generation (=0 is fixed by fBfix, not 0
  55                                    //impact parameter is generated in each event between fBfmin
  56                                    //and fBmax according with Glauber model (f-la 30)
  57   fHydjetParams.fBfix=0.;
  58   fHydjetParams.fBmin=0.;
  59   fHydjetParams.fBmax=0.5; //0-5% centrality
  60   fHydjetParams.fT = 0.170;
  61   fHydjetParams.fMuB = 0.0;
  62   fHydjetParams.fMuS = 0.0;
  63   fHydjetParams.fMuI3 = 0.0;
  64   fHydjetParams.fThFO = 0.130;
  65   fHydjetParams.fMu_th_pip = 0.0;
  66   fHydjetParams.fSeed=0;
  67   fHydjetParams.fTau=10.;
  68   fHydjetParams.fSigmaTau=3.;
  69   fHydjetParams.fR=11.;
  70   fHydjetParams.fYlmax=4.0;
  71   fHydjetParams.fUmax=1.1;
  72   fHydjetParams.fDelta=0.;
  73   fHydjetParams.fEpsilon=0.;
  74   fHydjetParams.fWeakDecay=0; //>=0 on ,-1 off
  75   fHydjetParams.fEtaType=1;//gaus
  76   fHydjetParams.fCorrS=1.;
  77   fHydjetParams.fNhsel=2;
  78   fHydjetParams.fIshad=1;
  79   fHydjetParams.fPtmin=7.0;
  80   fHydjetParams.fT0=0.8;
  81   fHydjetParams.fTau0=0.1;
  82   fHydjetParams.fNf=0;
  83   fHydjetParams.fIenglu=0;
  84   fHydjetParams.fIanglu=0;
  85
  86
  87   // Set reasonable default values for RHIC
  88 /*
  89   fHydjetParams.fSqrtS=200; //RHIC
  90   fHydjetParams.fAw=197;//Au-Au
  91   fHydjetParams.fIfb    = 1;
  92   fHydjetParams.fBfix=0.;
  93   fHydjetParams.fBmin=0.;
  94   fHydjetParams.fBmax=0.5; //0-5% centrality
  95   fHydjetParams.fT = 0.165;
  96   fHydjetParams.fMuB = 0.0285;
  97   fHydjetParams.fMuS = 0.007;
  98   fHydjetParams.fMuI3 = -0.001;
  99   fHydjetParams.fThFO = 0.100;
 100   fHydjetParams.fMu_th_pip = 0.053;
 101   fHydjetParams.fSeed=0;
 102   fHydjetParams.fTau=8.;
 103   fHydjetParams.fSigmaTau=2.;
 104   fHydjetParams.fR=10.;
 105   fHydjetParams.fYlmax=3.3;
 106   fHydjetParams.fUmax=1.1;
 107   fHydjetParams.fDelta=0.;
 108   fHydjetParams.fEpsilon=0.;
 109   fHydjetParams.fWeakDecay=0;//>=0 on ,-1 off
 110   fHydjetParams.fEtaType=1;//gaus
 111   fHydjetParams.fCorrS=1.;
 112   fHydjetParams.fNhsel=2;
 113   fHydjetParams.fIshad=1;
 114   fHydjetParams.fPtmin=3.4;
 115   fHydjetParams.fT0=0.3;
 116   fHydjetParams.fTau0=0.4;
 117   fHydjetParams.fNf=2;
 118   fHydjetParams.fIenglu=0;
 119   fHydjetParams.fIanglu=0;
 120 */
 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   //  cout << "TUHKMgen::TUHKMgen() OUT" << endl;
 130 }
 131
 132 //______________________________________________________________________________
 133 TUHKMgen::~TUHKMgen()
 134 {
 135   // destructor, deletes the InitialStateHydjet object
 136   if(fInitialState)
 137     delete fInitialState;
 138 }
 139
 140 void TUHKMgen::SetAllParametersRHIC()
 141 {
 142   // Set reasonable input parameters for 0-5% central Au+Au collisions
 143   // at 200 GeV at RHIC
 144   SetEcms(200.0);                // RHIC top energy
 145   SetAw(197);                    // Au+Au
 146   SetIfb(1);
 147   SetBfix(0.);
 148   SetBmin(0.0);                  // 0%
 149   SetBmax(0.5);                  // 5%
 150   SetChFrzTemperature(0.165);    // T_ch = 165 MeV
 151   SetMuB(0.0285);                // mu_B = 28.5 MeV
 152   SetMuS(0.007);                 // mu_S = 7 MeV
 153   SetMuQ(-0.001);                // mu_Q = -1 MeV
 154   SetThFrzTemperature(0.100);    // T_th = 100 MeV
 155   SetMuPionThermal(0.053);       // mu_th_pion = 53 MeV
 156   SetSeed(0);                    // use UNIX time
 157   SetTauB(8.0);                  // tau = 8 fm/c
 158   SetSigmaTau(2.0);              // sigma_tau = 2 fm/c
 159   SetRmaxB(10.0);                // fR = 10 fm
 160   SetYlMax(3.3);                 // fYmax = 3.3
 161   SetEtaRMax(1.1);               // Umax = 1.1
 162   SetMomAsymmPar(0.0);           // delta = 0.0
 163   SetCoordAsymmPar(0.0);         // epsilon = 0.0
 164   SetFlagWeakDecay(0.0);           // weak decay on (<0 off !!!)
 165   SetEtaType(1);                 // gaus distributed with fYmax dispersion (0 means boost invariant)
 166   SetGammaS(1.0);                // gammaS = 1.0 (no strangeness canonical suppresion)
 167   SetPyquenNhsel(2);             // hydro on, jets on, jet quenching on
 168   SetPyquenShad(1);              // shadowing on (0 off)
 169   SetPyquenPtmin(3.4);           // ptmin = 3.4 GeV/c
 170   SetPyquenT0(0.3);              // T0 = 300 MeV
 171   SetPyquenTau0(0.4);            // tau0 = 0.4 fm/c
 172   SetPyquenNf(2);                // 2 flavours
 173   SetPyquenIenglu(0);            // radiative and collisional energy loss
 174   SetPyquenIanglu(0);            // small gluon angular distribution
 175 }
 176
 177 void TUHKMgen::SetAllParametersLHC()
 178 {
 179   // Set reasonable input parameters for 0-5% Pb+Pb collisions at 5.5 TeV at LHC
 180   SetEcms(5500.0);               // LHC
 181   SetAw(207);                    // Pb+Pb
 182   SetIfb(1);
 183   SetBfix(0.);                  // 0
 184   SetBmin(0.0);                  // 0%
 185   SetBmax(0.5);                  // 5%
 186   SetChFrzTemperature(0.170);    // T_ch = 170 MeV
 187   SetMuB(0.0);                   // mu_B = 0 MeV
 188   SetMuS(0.0);                   // mu_S = 0 MeV
 189   SetMuQ(0.0);                   // mu_Q = 0 MeV
 190   SetThFrzTemperature(0.130);    // T_th = 130 MeV
 191   SetMuPionThermal(0.0);         // mu_th_pion = 0 MeV
 192   SetSeed(0);                    // use UNIX time
 193   SetTauB(10.0);                 // tau = 10 fm/c
 194   SetSigmaTau(3.0);              // sigma_tau = 3 fm/c
 195   SetRmaxB(11.0);                // fR = 11 fm
 196   SetYlMax(4.0);                 // fYmax = 4.0
 197   SetEtaRMax(1.1);               // Umax = 1.1
 198   SetMomAsymmPar(0.0);           // delta = 0.0
 199   SetCoordAsymmPar(0.0);         // epsilon = 0.0
 200   SetFlagWeakDecay(0);           // weak decay on (<0 off !!!)
 201   SetEtaType(1);                 // gaus distributed with fYmax dispersion (0 means boost invariant)
 202   SetGammaS(1.0);                // gammaS = 1.0 (no strangeness canonical suppresion)
 203   SetPyquenNhsel(2);             // hydro on, jets on, jet quenching on
 204   SetPyquenShad(1);              // shadowing on (0 off)
 205   SetPyquenPtmin(7.0);           // ptmin = 7.0 GeV/c
 206   SetPyquenT0(0.8);              // T0 = 800 MeV
 207   SetPyquenTau0(0.1);            // tau0 = 0.4 fm/c
 208   SetPyquenNf(0);                // 0 flavours
 209   SetPyquenIenglu(0);            // radiative and collisional energy loss
 210   SetPyquenIanglu(0);            // small gluon angular distribution
 211 }
 212
 213 TObjArray* TUHKMgen::ImportParticles(const Option_t *)
 214 {
 215   // Function overloading the TGenerator::ImportParticles() member function.
 216   // The particles from the local particle list (fSecondariesList) are
 217   // forwarded to the TGenerator::fParticles
 218   fParticles->Clear();
 219   Int_t nump = 0;
 220   LPIT_t it,e;
 221
 222   for(it = fSecondariesList.begin(), e = fSecondariesList.end(); it != e; it++) {
 223     TVector3 pos(it->Pos().Vect());
 224     TVector3 mom(it->Mom().Vect());
 225     Float_t m1 = it->TableMass();
 226     Int_t im1 = it->GetMother();
 227     Int_t im2 = -1;
 228     Int_t id1 = -1;
 229     Int_t id2 = -1;
 230
 231     Int_t type = it->GetType();  // 0-hydro, 1-jets
 232
 233     if (im1> -1) {
 234       TParticle *mother = (TParticle*) (fParticles->UncheckedAt(im1+1));
 235       mother->SetLastDaughter(nump);
 236       if (mother->GetFirstDaughter()==-1)
 237         mother->SetFirstDaughter(nump+1);
 238     }
 239
 240     nump++;
 241
 242     TParticle* p = new TParticle(it->Encoding(), type,                          //pdg,stat
 243                                  im1, im2, id1, id2,                                            //m1,m2,d1,d2
 244                                  mom[0], mom[1], mom[2], TMath::Sqrt(mom.Mag2() + m1 * m1),     //px,py,pz,e
 245                                  pos[0]*1.e-13, pos[1]*1.e-13, pos[2]*1.e-13,
 246                                  it->T()*1.e-13/3e10);                          //x,y,z,t
 247
 248      p->SetUniqueID(nump);
 249      fParticles->Add(p);
 250   } //end for
 251
 252   fAllocator.FreeList(fSecondariesList);
 253
 254   return fParticles;
 255 }
 256
 257 Int_t TUHKMgen::ImportParticles(TClonesArray *particles, const Option_t* option)
 258 {
 259   // Function overloading the TGenerator::ImportParticles() member function.
 260   // The particles from the local particle list (fSecondariesList) are
 261   // forwarded to the TGenerator::fParticles
 262
 263
 264
 265   if(particles==0) return 0;
 266   TClonesArray &particlesR=*particles;
 267   particlesR.Clear();
 268
 269   Int_t numprim,numsec;  numprim=numsec=0;
 270   Int_t nump = 0;
 271   LPIT_t it,e;
 272
 273   for(it = fSecondariesList.begin(), e = fSecondariesList.end(); it != e; it++) {
 274     TVector3 pos(it->Pos().Vect());
 275     TVector3 mom(it->Mom().Vect());
 276     Float_t m1 = it->TableMass();
 277     Int_t im1 = it->GetMother();
 278     Int_t im2 = -1;
 279     Int_t id1 = -1;
 280     Int_t id2 = -1;
 281
 282     Int_t type = it->GetType();
 283
 284     if (im1> -1) {
 285      // particle not a primary -> set the daughter indexes for the mother particle"<< endl;
 286       TParticle *mother = (TParticle*) (particlesR.UncheckedAt(im1));
 287       mother->SetLastDaughter(nump);
 288       if(mother->GetFirstDaughter()==-1)
 289         mother->SetFirstDaughter(nump);
 290     }
 291
 292
 293     new (particlesR[nump]) TParticle(it->Encoding(), type,                                              //pdg,stat
 294                                      im1, im2, id1, id2,                                                //m1,m2,d1,d2
 295                                      mom[0], mom[1], mom[2], TMath::Sqrt(mom.Mag2() + m1 * m1), //px,py,pz,e
 296                                      pos[0]*1.e-13, pos[1]*1.e-13, pos[2]*1.e-13,
 297                                      it->T()*1.e-13/3e10);                      //x,y,z,t
 298
 299     particlesR[nump]->SetUniqueID(nump);
 300     nump++;
 301     numsec++;
 302   }//end for
 303
 304   fSecondariesList.clear();
 305   printf("Scan and add prim %d sec %d and all %d particles\n",
 306          numprim,numsec,nump);
 307   return nump;
 308 }
 309
 310 //______________________________________________________________________________
 311 void TUHKMgen::Initialize()
 312 {
 313   // Function overloading the TGenerator::Initialize() member function.
 314   // The Monte-Carlo model is initialized (input parameters are transmitted,
 315   // particle list and decay channels are loaded, average multiplicities are calculated, etc.)
 316
 317   fInitialState = new InitialStateHydjet();
 318   SetAllParameters();
 319   fInitialState->LoadPDGInfo();
 320   // set the stable flags
 321   for(Int_t i=0; i<fStableFlagged; i++)
 322     fInitialState->SetPDGParticleStable(fStableFlagPDG[i], fStableFlagStatus[i]);
 323
 324   if(!fInitialState->MultIni())
 325     Error("TUHKMgen::Initialize()", "Bad status return from MultIni(). Check it out!! \n");  //
 326 }
 327
 328 void TUHKMgen::Print(const Option_t*) const
 329 {
 330   cout << "TUHKMgen::Print() method not implemented yet!!" << endl;
 331 }
 332
 333 void TUHKMgen::GenerateEvent()
 334 {
 335   // Member function overloading the TGenerator::GenerateEvent()
 336   // The HYDJET++ model is run and the particle lists (fSourceList and fSecondariesList) are filled
 337
 338   fInitialState->Initialize(fSecondariesList, fAllocator);
 339
 340   if(fSecondariesList.empty())Error("TUHKMgen::GenerateEvent()", "Source particle list empty after fireball initialization!! \n");  //
 341
 342   // Run the decays
 343
 344   if(fInitialState->RunDecays())
 345     fInitialState->Evolve(fSecondariesList, fAllocator, fInitialState->GetWeakDecayLimit());
 346
 347 }
 348
 349 void TUHKMgen::SetAllParameters() {
 350   // forward all model input parameters to the InitialStateHydjet object
 351   // which will handle all the Monte-Carlo simulation using HYDJET++ model
 352
 353   fInitialState->fParams.fSqrtS = fHydjetParams.fSqrtS;
 354   fInitialState->fParams.fAw = fHydjetParams.fAw;
 355   fInitialState->fParams.fIfb = fHydjetParams.fIfb ;
 356   fInitialState->fParams.fBfix = fHydjetParams.fBfix;
 357   fInitialState->fParams.fBmin = fHydjetParams.fBmin;
 358   fInitialState->fParams.fBmax = fHydjetParams.fBmax;
 359   fInitialState->fParams.fSeed = fHydjetParams.fSeed;
 360   fInitialState->fParams.fT = fHydjetParams.fT;
 361   fInitialState->fParams.fMuB = fHydjetParams.fMuB;
 362   fInitialState->fParams.fMuS = fHydjetParams.fMuS;
 363   fInitialState->fParams.fMuI3 = fHydjetParams.fMuI3;
 364   fInitialState->fParams.fThFO = fHydjetParams.fThFO;
 365   fInitialState->fParams.fMu_th_pip = fHydjetParams.fMu_th_pip;
 366
 367   fInitialState->fParams.fTau = fHydjetParams.fTau;
 368   fInitialState->fParams.fSigmaTau = fHydjetParams.fSigmaTau;
 369   fInitialState->fParams.fR = fHydjetParams.fR;
 370   fInitialState->fParams.fYlmax = fHydjetParams.fYlmax;
 371   fInitialState->fParams.fUmax = fHydjetParams.fUmax;
 372   fInitialState->fParams.fDelta = fHydjetParams.fDelta;
 373   fInitialState->fParams.fEpsilon = fHydjetParams.fEpsilon;
 374
 375   fInitialState->fParams.fWeakDecay = fHydjetParams.fWeakDecay;
 376   fInitialState->fParams.fDecay = 1;                   // run decays
 377   fInitialState->fParams.fCorrS = fHydjetParams.fCorrS;
 378   fInitialState->fParams.fEtaType = fHydjetParams.fEtaType;
 379
 380   fInitialState->fParams.fPtmin = fHydjetParams.fPtmin;
 381   fInitialState->fParams.fNhsel = fHydjetParams.fNhsel;
 382   fInitialState->fParams.fIshad = fHydjetParams.fIshad;
 383   fInitialState->fParams.fT0 = fHydjetParams.fT0;
 384   fInitialState->fParams.fTau0 = fHydjetParams.fTau0;
 385   fInitialState->fParams.fNf = fHydjetParams.fNf;
 386   fInitialState->fParams.fIenglu = fHydjetParams.fIenglu;
 387   fInitialState->fParams.fIanglu = fHydjetParams.fIanglu;
 388
 389   fInitialState->SetPDGParticleFilename(fParticleFilename);
 390   fInitialState->SetPDGDecayFilename(fDecayFilename);
 391   //  fInitialState->SetUseCharmParticles(fUseCharmParticles);
 392   //  fInitialState->SetWidthRange(fMinWidth, fMaxWidth);
 393   //  fInitialState->SetMassRange(fMinMass, fMaxMass);
 394   //  for(Int_t i=0; i<fStableFlagged; i++)
 395   //    fInitialState->SetPDGParticleStable(fStableFlagPDG[i], fStableFlagStatus[i]);
 396   //  cout << "TUHKMgen::SetAllParameters() OUT" << endl;
 397 }
 398