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