Corrected warnings, coding conventions, bug fixes (Ionut)
authorhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 26 Jun 2009 12:48:29 +0000 (12:48 +0000)
committerhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 26 Jun 2009 12:48:29 +0000 (12:48 +0000)
19 files changed:
TUHKMgen/AliGenUHKM.cxx
TUHKMgen/AliGenUHKM.h
TUHKMgen/PYQUEN/progs_fortran.f
TUHKMgen/TUHKMgen.cxx
TUHKMgen/TUHKMgen.h
TUHKMgen/UHKM/DatabasePDG.cxx
TUHKMgen/UHKM/DatabasePDG.h
TUHKMgen/UHKM/GrandCanonical.h
TUHKMgen/UHKM/HadronDecayer.cxx
TUHKMgen/UHKM/HadronDecayer.h
TUHKMgen/UHKM/InitialState.cxx
TUHKMgen/UHKM/InitialState.h
TUHKMgen/UHKM/InitialStateHydjet.cxx
TUHKMgen/UHKM/InitialStateHydjet.h
TUHKMgen/UHKM/Particle.cxx
TUHKMgen/UHKM/Particle.h
TUHKMgen/UHKM/ParticlePDG.h
TUHKMgen/UHKM/StrangeDensity.h
TUHKMgen/UHKM/StrangePotential.h

index 42a89ba..d03f251 100755 (executable)
 #include <string>
 using namespace std;
 
-//class TClonesArray;
-//enum TMCProcess;
-//class AliRun;
-//class TSystem;
-//class AliDecayer; 
 
 ClassImp(AliGenUHKM)
 
@@ -59,7 +54,6 @@ AliGenUHKM::AliGenUHKM()
 {
   // Default constructor setting up default reasonable parameter values
   // for central Pb+Pb collisions at 5.5TeV 
-  //  cout << "AliGenUHKM::AliGenUHKM() IN" << endl;
   
   // LHC
   fHydjetParams.fSqrtS=5500; //LHC
@@ -132,7 +126,6 @@ AliGenUHKM::AliGenUHKM()
   }
   fStableFlagged = 0;
 
-  //  cout << "AliGenUHKM::AliGenUHKM() OUT" << endl;
 }
 
 //_______________________________________
@@ -146,7 +139,7 @@ AliGenUHKM::AliGenUHKM(Int_t npart)
   // Constructor specifying the size of the particle table
   // and setting up default reasonable parameter values
   // for central Pb+Pb collisions at 5.5TeV
-  //  cout << "AliGenUHKM::AliGenUHKM(Int_t) IN" << endl;
+
   fName = "UHKM";
   fTitle= "Particle Generator using UHKM 3.0";
 
@@ -223,13 +216,13 @@ AliGenUHKM::AliGenUHKM(Int_t npart)
   }
   fStableFlagged = 0;  
 
-  //  cout << "AliGenUHKM::AliGenUHKM(Int_t) OUT" << endl;
 }
 
 //__________________________________________
 AliGenUHKM::~AliGenUHKM()
 {
   // Destructor, do nothing
+  //  delete fParticles;
 }
 
 void AliGenUHKM::SetAllParametersRHIC() 
@@ -310,19 +303,15 @@ void AliGenUHKM::Init()
   // further to the model.
   // HYDJET++ is initialized (average multiplicities are calculated, particle species definitions and decay
   // channels are loaded, etc.)
-  //  cout << "AliGenUHKM::Init() IN" << endl;
 
   SetMC(new TUHKMgen());
   fUHKMgen = (TUHKMgen*) fMCEvGen;
   SetAllParameters();
   
   AliGenMC::Init();
-
   fUHKMgen->Initialize();
-  CheckPDGTable();    //
-
-  fUHKMgen->Print();  //
-  //  cout << "AliGenUHKM::Init() OUT" << endl;
+  CheckPDGTable();    
+  fUHKMgen->Print();  
 }
 
 //________________________________________
@@ -330,7 +319,7 @@ void AliGenUHKM::Generate()
 {
   // Generate one HYDJET++ event, get the output and push particles further 
   // to AliRoot's stack 
-  //  cout << "AliGenUHKM::Generate() IN" << endl;
+
   Float_t polar[3] = {0,0,0};
   Float_t origin[3]   = {0,0,0};
   Float_t origin0[3]  = {0,0,0};
@@ -341,7 +330,6 @@ void AliGenUHKM::Generate()
   for(Int_t j=0; j<3; j++) origin0[j] = fVertex[j];
 
   fTrials = 0;
- // cout << "AliGenUHKM::Generate() fTrials = " << fTrials << endl;
 
   Int_t nt  = 0;
 
@@ -351,7 +339,6 @@ void AliGenUHKM::Generate()
   fUHKMgen->ImportParticles(&fParticles,"All");
 
   Int_t np = fParticles.GetEntriesFast();
-  //  cout << "AliGenUHKM::Generate() GetEntries  " <<np<< endl;
 
 
   Int_t* idsOnStack = new Int_t[np];
@@ -359,33 +346,17 @@ void AliGenUHKM::Generate()
   for(Int_t i=0; i<np; i++) newPos[i] = i;
 
   //_________ Loop for particle selection
-  //  for(Int_t i=1; i<np; i++) {
-  for(Int_t i=1; i<np; i++) {
+  for(Int_t i=0; i<np; i++) {
     // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     // the particle indexes are 0 based but fParticles is a 1 based array
     // -1 is the trivial code (when it does not exist)
     // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     TParticle *iparticle = (TParticle*)fParticles.At(i);
-    //   cout << "AliGenUHKM::Generate() particle #" << i << " in fParticles *********************"<< endl;
-
     Int_t kf = iparticle->GetPdgCode();
-    //    if(kf==130 || kf==-130)
-    //      cout << "AliGenUHKM::Generate() PDG = " << kf << endl;
-
     Bool_t hasMother = (iparticle->GetFirstMother() >= 0);
-
-    //    cout << "AliGenUHKM::Generate() mother index in fParticles = " 
-    //  << (iparticle->GetFirstMother()==-1 ? -1 : iparticle->GetFirstMother()+1)  
-    //  << " ; hasMother = " << hasMother << endl;
-    //    cout << "AliGenUHKM::Generate() type (0-hydro;1-jet) = " << iparticle->GetStatusCode() << endl;
     Bool_t hasDaughter = (iparticle->GetNDaughters() > 0);
 
-    // cout << "AliGenUHKM::Generate() n.daughters = " << iparticle->GetNDaughters() 
-    //<< " ; hasDaughter = " << hasDaughter << endl;
-
-
     if(hasDaughter) {
-      //      cout << "AliGenUHKM::Generate() decayed particle (not trackable)" << endl;
       // This particle has decayed
       // It will not be tracked
       // Add it only once with coordinates not
@@ -400,28 +371,22 @@ void AliGenUHKM::Generate()
       v[2] = iparticle->Vz();
       Float_t time = iparticle->T();
 
-      //      if(kf==130 || kf==-130) cout << "type = " << iparticle->GetStatusCode() << endl; //j1/h0
-
       Int_t imo = -1;
       if(hasMother) {
         imo = iparticle->GetFirstMother(); //index of mother particle in fParticles
       } // if has mother
       Bool_t trackFlag = kFALSE;   // tFlag is kFALSE --> do not track the particle
 
-      //      printf("Pushing Track %d with status %d mother %d\n", kf, trackFlag, imo>=0?idsOnStack[imo]:imo);
       PushTrack(trackFlag, (imo>=0 ? idsOnStack[imo+1] : imo), kf,
                 p[0], p[1], p[2], energy,
                 v[0], v[1], v[2], time,
                 polar[0], polar[1], polar[2],
                 (hasMother ? kPDecay : kPNoProcess), nt);
-      //      cout << "AliGenUHKM::Generate() pushed on stack with stack index = " << nt 
-      //          << "; mother index on stack = " << (imo>=0 ? idsOnStack[imo+1] : imo) << endl;
       idsOnStack[i] = nt;
       fNprimaries++;
       KeepTrack(nt);
     }
     else {
-      //            cout << "AliGenUHKM::Generate() final particle --> push it twice on the stack" << endl;
       // This is a final state particle
       // It will be tracked
       // Add it TWICE to the stack !!!
@@ -439,7 +404,6 @@ void AliGenUHKM::Generate()
       v[2] = iparticle->Vz();
 
       Int_t type    = iparticle->GetStatusCode(); // 1-from jet / 0-from hydro 
-      //if(kf==130 || kf==-130) cout << "type = " << iparticle->GetStatusCode() << endl; //j1/h0
       Int_t coeffT=1;
       if(type==1) coeffT=-1; //to separate particles from jets
 
@@ -449,13 +413,11 @@ void AliGenUHKM::Generate()
         imo = iparticle->GetFirstMother();
       } // if has mother
       Bool_t trackFlag = kFALSE;  // tFlag = kFALSE --> do not track this one, its for femtoscopy
-      PushTrack(trackFlag, (imo>=0 ? idsOnStack[imo+1] : imo), kf,
+      PushTrack(trackFlag, (imo>=0 ? idsOnStack[imo] : imo), kf,
                 p[0], p[1], p[2], energy,
                 v[0], v[1], v[2], (iparticle->T())*coeffT,
                 polar[0], polar[1], polar[2],
                 hasMother ? kPDecay:kPNoProcess, nt);
-      //      cout << "AliGenUHKM::Generate() pushed on stack with stack index = " << nt
-      //           << "; mother index on stack = " << (imo>=0 ? idsOnStack[imo+1] : imo) << endl;
 
       idsOnStack[i] = nt;
       fNprimaries++;
@@ -467,15 +429,12 @@ void AliGenUHKM::Generate()
       imo = nt;
       
       trackFlag = kTRUE;    // tFlag = kTRUE --> track this one
-      //      cout << "AliGenUHKM::Generate() trackFlag = " << trackFlag << endl;
 
       PushTrack(trackFlag, imo, kf,
                 p[0], p[1], p[2], energy,
                 origin[0], origin[1], origin[2], iparticle->T(),
                 polar[0], polar[1], polar[2],
                 hasMother ? kPDecay:kPNoProcess, nt);
-      //      cout << "AliGenUHKM::Generate() pushed on stack with stack index = " << nt
-      //           << "; mother index on stack = " << imo << endl;
       fNprimaries++;
       KeepTrack(nt);
     }
@@ -508,10 +467,6 @@ void AliGenUHKM::Generate()
 
   delete idsOnStack;
 
-  //  gAlice->SetGenEventHeader(header);
-
-  //  printf(" Finish Generate .. %d ..\n",nt);
-  //  cout << "AliGenUHKM::Generate() OUT" << endl;
 }
 
 void AliGenUHKM::Copy(TObject &) const
@@ -522,8 +477,6 @@ void AliGenUHKM::Copy(TObject &) const
 void AliGenUHKM::SetAllParameters() {
   // Forward all input parameters to the TGenerator::TUHKMgen object
 
-  //  cout << "AliGenUHKM::SetAllParameters() IN" << endl;
-
   fUHKMgen->SetEcms(fHydjetParams.fSqrtS);
   fUHKMgen->SetBmin(fHydjetParams.fBmin);
   fUHKMgen->SetBmax(fHydjetParams.fBmax);
@@ -546,7 +499,6 @@ void AliGenUHKM::SetAllParameters() {
   fUHKMgen->SetCoordAsymmPar(fHydjetParams.fEpsilon);
 
   fUHKMgen->SetGammaS(fHydjetParams.fCorrS);
-  // fUHKMgen->SetHdec(fHydjetParams.fTime);
   fUHKMgen->SetEtaType(fHydjetParams.fEtaType);
   fUHKMgen->SetFlagWeakDecay(fHydjetParams.fWeakDecay);
 
@@ -606,7 +558,6 @@ cout<<"-----Volume parameters -------------- "<<endl;
 cout<<" --------Flags------ "<<endl;
 
  cout<<" GammaS= "<<fHydjetParams.fCorrS<<endl;
-  // fUHKMgen->SetHdec(fHydjetParams.fTime<<endl;
  cout<<" EtaType= "<<fHydjetParams.fEtaType<<endl;
  cout<<" FlagWeakDecay= "<<fHydjetParams.fWeakDecay<<endl;
 
@@ -637,19 +588,14 @@ void AliGenUHKM::CheckPDGTable() {
   // Add temporarely all particle definitions from HYDJET++ which miss in the ROOT's PDG tables
   // to the TDatabasePDG table.
 
-  //  cout << "AliGenUHKM::CheckPDGTable()   IN" << endl;
-  //TDabasePDG *rootPDG  = TDatabasePDG::Instance();         // ROOT's PDG table
   DatabasePDG *uhkmPDG = fUHKMgen->PDGInfo();              // UHKM's PDG table
   TParticlePDG *rootTestParticle;
   ParticlePDG *uhkmTestParticle;
   
-  //  cout << "particles with good status in UHKM table = " << uhkmPDG->GetNParticles() << endl;
   // loop over all particles in the SHARE table
   for(Int_t i=0; i<uhkmPDG->GetNParticles(); i++) {
-    //    cout << "particle #" << i << " ================" << endl;
     // get a particle specie
     uhkmTestParticle = uhkmPDG->GetPDGParticleByIndex(i);
-    //    cout << "PDG = " << uhkmTestParticle->GetPDG() << endl;
     // check if this code exists in ROOT's table
     rootTestParticle = TDatabasePDG::Instance()->GetParticle(uhkmTestParticle->GetPDG());
     if(!rootTestParticle) {    // if not then add it to the ROOT's PDG database
@@ -660,14 +606,10 @@ void AliGenUHKM::CheckPDGTable() {
                                            uhkmTestParticle->GetWidth(), 
                                            (Int_t(uhkmTestParticle->GetBaryonNumber())==0 ? "meson" : "baryon"),
                                            uhkmTestParticle->GetPDG());    
-      //      cout << "Not found in ROOT's PDG database --> added now" << endl;
       if(uhkmTestParticle->GetWidth()<1e-10) 
        cout << uhkmTestParticle->GetPDG() << "its mass "
             << TDatabasePDG::Instance()->GetParticle(uhkmTestParticle->GetPDG())->Mass()
             << TDatabasePDG::Instance()->GetParticle(uhkmTestParticle->GetPDG())->Width() << endl;  
     }
-    //    else
-      //      cout << "Found in ROOT's PDG database --> not added" << endl;
   }  // end for
-  //  cout << "AliGenUHKM::CheckPDGTable()   OUT" << endl;
 }
index 07e3cf9..546a8ce 100755 (executable)
@@ -9,7 +9,7 @@
 #include "AliGenMC.h"
 #include <TString.h>
 #include "TUHKMgen.h"
-#ifndef INITIALSTATEHYDJET_INCLUDED
+#ifndef INITIALSTATEHYDJET_H
 #include "InitialStateHydjet.h"
 #endif
 #include "TParticle.h"
@@ -17,9 +17,6 @@
 #include <string>
 using namespace std;
 
-//class TUHKMgen;
-//class TParticle;
-//class InitialStateHydjet;
 
 class AliGenUHKM : public AliGenMC
 {
index e172238..cab5612 100644 (file)
@@ -316,10 +316,7 @@ c-ml
 
 c--      save/lujets/,/hyjets/,/hyipar/,/hyfpar/,/hyflow/,/hyjpar/,/hypyin/
 
-
-* reset lujets and hyjets arrays before event generation 
-
-         write(*,*)'in hyevnt 0'
+* reset hyjets arrays before event generation 
 
       n=0 
       nhj=0 
@@ -410,7 +407,7 @@ c       end do
 c      end do
       
 c-ml
-c-      call luedit(2)                      ! remove unstable particles and partons 
+c      call luedit(2)                      ! remove unstable particles and partons 
 
 
 
index 47361a8..2796f9a 100755 (executable)
@@ -37,22 +37,15 @@ ClassImp(TUHKMgen)
 TUHKMgen::TUHKMgen() : 
   TGenerator("UHKM","UHKM"),
   fInitialState(0x0),
-  fAllocator(),
-  fSourceList(),
-  fSecondariesList(),
   fNPprim(0),
   fNPsec(0),
   fHydjetParams(),
   fStableFlagged(0)
-  //  fUseCharmParticles(kFALSE),
-  //  fMinWidth(0.0),
-  //  fMaxWidth(10.0),
-  //  fMinMass(0.0001),
-  //  fMaxMass(10.0)
 {
   // default constructor setting reasonable defaults for initial parameters (central Pb+Pb at 5.5 TeV)
 
-  //  cout << "TUHKMgen::TUHKMgen() IN" << endl;
+  ParticleAllocator fAllocator;
+  List_t fSecondariesList;
 
   // Set reasonable default values for LHC
   
@@ -140,10 +133,8 @@ TUHKMgen::TUHKMgen() :
 TUHKMgen::~TUHKMgen()
 {
   // destructor, deletes the InitialStateHydjet object
-  //  cout << "TUHKMgen::~TUHKMgen() IN" << endl;
   if(fInitialState)
     delete fInitialState;
-  //  cout << "TUHKMgen::~TUHKMgen() OUT" << endl;
 }
 
 void TUHKMgen::SetAllParametersRHIC()
@@ -224,13 +215,11 @@ TObjArray* TUHKMgen::ImportParticles(const Option_t *)
   // Function overloading the TGenerator::ImportParticles() member function.
   // The particles from the local particle list (fSecondariesList) are
   // forwarded to the TGenerator::fParticles 
-  //  cout << "TObjArray* TUHKMgen::ImportParticles(Option_t) IN" << endl;
   fParticles->Clear();
-  //  cout << "TObjArray* TUHKMgen::ImportParticles(Option_t): UHKM stack contains " << fNPsec << " particles." << endl;
   Int_t nump = 0;
   LPIT_t it,e;
    
-  for(it = fSecondariesList.begin(), e = fSecondariesList.end(); it != e; ++it) {
+  for(it = fSecondariesList.begin(), e = fSecondariesList.end(); it != e; it++) {
     TVector3 pos(it->Pos().Vect());
     TVector3 mom(it->Mom().Vect());
     Float_t m1 = it->TableMass();
@@ -239,15 +228,13 @@ TObjArray* TUHKMgen::ImportParticles(const Option_t *)
     Int_t id1 = -1;
     Int_t id2 = -1;
 
-    //    Int_t nd = it->GetNDaughters();
-    
     Int_t type = it->GetType();  // 0-hydro, 1-jets
 
     if (im1> -1) {
-      TParticle *mother = (TParticle*) (fParticles->UncheckedAt(im1));    
+      TParticle *mother = (TParticle*) (fParticles->UncheckedAt(im1+1));          
       mother->SetLastDaughter(nump);
       if (mother->GetFirstDaughter()==-1)
-       mother->SetFirstDaughter(nump);
+       mother->SetFirstDaughter(nump+1);
     }
 
     nump++;
@@ -262,9 +249,8 @@ TObjArray* TUHKMgen::ImportParticles(const Option_t *)
      fParticles->Add(p);
   } //end for
  
-  fAllocator.FreeList(fSourceList);
   fAllocator.FreeList(fSecondariesList);
-  //  cout << "TObjArray* TUHKMgen::ImportParticles(Option_t) OUT" << endl;
+
   return fParticles;
 }
 
@@ -273,69 +259,51 @@ Int_t TUHKMgen::ImportParticles(TClonesArray *particles, const Option_t* option)
   // Function overloading the TGenerator::ImportParticles() member function.
   // The particles from the local particle list (fSecondariesList) are
   // forwarded to the TGenerator::fParticles
-  //  cout << "Int_t TUHKMgen::ImportParticles(TClonesArray*, Option_t) IN" << endl;
+
+  
+
   if(particles==0) return 0;
   TClonesArray &particlesR=*particles;
   particlesR.Clear();
-
  
   Int_t numprim,numsec;  numprim=numsec=0;
   Int_t nump = 0;
   LPIT_t it,e;
   
-    cout << "TUHKMgen::ImportParticles() option(All or Sec) = " << option << endl;
-  for(it = fSecondariesList.begin(), e = fSecondariesList.end(); it != e; ++it) {
-    //!!! for(Int_t pp=0;pp<100;pp++) {
-       //      cout << "TUHKMgen::ImportParticles() import particle pdg(" << it->Encoding() << ")" << endl;
-    TVector3 pos(it->Pos().Vect());  //
-    //!!! TVector3 pos(0.0, 0.1, 0.2);
-    TVector3 mom(it->Mom().Vect());  //
-    //!!! TVector3 mom(-0.1, 0.1, 0.2);
-    Float_t m1 = it->TableMass();  //
-    //!!!    Float_t m1 = 0.139;
-    Int_t im1 = it->GetMother();  //
-    //!!!    Int_t im1 = (pp<50? -1 : pp-50);
-    //    Int_t nd = 0;
+  for(it = fSecondariesList.begin(), e = fSecondariesList.end(); it != e; it++) {
+    TVector3 pos(it->Pos().Vect());  
+    TVector3 mom(it->Mom().Vect());  
+    Float_t m1 = it->TableMass();  
+    Int_t im1 = it->GetMother();  
     Int_t im2 = -1;
     Int_t id1 = -1;
     Int_t id2 = -1;
     
-    //    nd = it->GetNDaughters();
-
-    Int_t type = it->GetType();  //  // 0-hydro, 1-jets
-    //!!!    Int_t type = 0;
-    
+    Int_t type = it->GetType();  
+      
     if (im1> -1) {
-      TParticle *mother = (TParticle*) (particlesR.UncheckedAt(im1+1));           
+     // particle not a primary -> set the daughter indexes for the mother particle"<< endl;
+      TParticle *mother = (TParticle*) (particlesR.UncheckedAt(im1));
       mother->SetLastDaughter(nump);
-      if (mother->GetFirstDaughter()==-1)
-       mother->SetFirstDaughter(nump);
+      if(mother->GetFirstDaughter()==-1)
+       mother->SetFirstDaughter(nump);
     }
       
-    nump++;
-
+    
     new (particlesR[nump]) TParticle(it->Encoding(), type,                                             //pdg,stat
                                     im1, im2, id1, id2,                                                //m1,m2,d1,d2
                                     mom[0], mom[1], mom[2], TMath::Sqrt(mom.Mag2() + m1 * m1), //px,py,pz,e
                                     pos[0]*1.e-13, pos[1]*1.e-13, pos[2]*1.e-13, 
                                     it->T()*1.e-13/3e10);                      //x,y,z,t
     
-    /*!!! new (particlesR[nump]) TParticle(211, type,                                         //pdg,stat
-                                     im1, im2, id1, id2,                                                //m1,m2,d1,d2
-                                     mom[0], mom[1], mom[2], TMath::Sqrt(mom.Mag2() + m1 * m1), //px,py,pz,e
-                                     pos[0]*1.e-13, pos[1]*1.e-13, pos[2]*1.e-13,
-                                     0.0*1.e-13/3e10);                       //x,y,z,t
-    !!!*/
-
     particlesR[nump]->SetUniqueID(nump);
+    nump++;
     numsec++;
   }//end for
   
-  fAllocator.FreeList(fSourceList);  //
-  fAllocator.FreeList(fSecondariesList);  //
-  //  printf("Scan and add prim %d sec %d and all %d particles\n",
-  //    numprim,numsec,nump);
-  //  cout << "Int_t TUHKMgen::ImportParticles(TClonesArray*, Option_t) OUT" << endl;
+  fSecondariesList.clear();
+  printf("Scan and add prim %d sec %d and all %d particles\n",
+        numprim,numsec,nump);
   return nump;
 }
 
@@ -345,17 +313,16 @@ void TUHKMgen::Initialize()
   // Function overloading the TGenerator::Initialize() member function.
   // The Monte-Carlo model is initialized (input parameters are transmitted, 
   // particle list and decay channels are loaded, average multiplicities are calculated, etc.)
-  //  cout << "TUHKMgen::Initialize() IN" << endl;
-  fInitialState = new InitialStateHydjet();  //
-  SetAllParameters();  //
-  fInitialState->LoadPDGInfo();  //
+  fInitialState = new InitialStateHydjet();  
+  SetAllParameters();  
+  fInitialState->LoadPDGInfo();  
   // set the stable flags
-  for(Int_t i=0; i<fStableFlagged; i++) //
-    fInitialState->SetPDGParticleStable(fStableFlagPDG[i], fStableFlagStatus[i]); //
+  for(Int_t i=0; i<fStableFlagged; i++) 
+    fInitialState->SetPDGParticleStable(fStableFlagPDG[i], fStableFlagStatus[i]); 
 
-  if(!fInitialState->MultIni()) //
+  if(!fInitialState->MultIni()) 
     Error("TUHKMgen::Initialize()", "Bad status return from MultIni(). Check it out!! \n");  //
-  //  cout << "TUHKMgen::Initialize() OUT" << endl;
 }
 
 void TUHKMgen::Print(const Option_t*) const
@@ -367,32 +334,22 @@ void TUHKMgen::GenerateEvent()
 {
   // Member function overloading the TGenerator::GenerateEvent()
   // The HYDJET++ model is run and the particle lists (fSourceList and fSecondariesList) are filled
-  //  cout << "TUHKMgen::GenerateEvent() IN" << endl;
- // cout << "TUHKMgen::GenerateEvent() IN fSourceList " <<endl;
-  //  cout<< " fInitialStrate "<<fInitialState<<endl;  //
-  //" fSourceList "<< fSourceList<<" fAllocator "<< fAllocator <<endl;
   
-  // Generate the initial particles
-  fInitialState->Initialize(fSourceList, fAllocator);  //
-  //  cout << "TUHKMgen::fInitialState->Initialize" << endl;
+  fInitialState->Initialize(fSecondariesList, fAllocator);  
  
-  if(fSourceList.empty()) //
-    Error("TUHKMgen::GenerateEvent()", "Source particle list empty after fireball initialization!! \n");  //
+  if(fSecondariesList.empty())Error("TUHKMgen::GenerateEvent()", "Source particle list empty after fireball initialization!! \n");  //
 
   // Run the decays
 
-  //  cout << "after Initilize: get_time, weak_decay_limit" <<fInitialState->GetTime()<<fInitialState->GetWeakDecayLimit()<< endl;   //
-  if(fInitialState->GetTime() >= 0.)  //
-    fInitialState->Evolve(fSourceList, fSecondariesList, fAllocator, fInitialState->GetWeakDecayLimit());  //
-  //  cout << "TUHKMgen::GenerateEvent() OUT" << endl;
+  if(fInitialState->RunDecays())  
+    fInitialState->Evolve(fSecondariesList, fAllocator, fInitialState->GetWeakDecayLimit());  
+  
 }
 
 void TUHKMgen::SetAllParameters() {
   // forward all model input parameters to the InitialStateHydjet object
   // which will handle all the Monte-Carlo simulation using HYDJET++ model
 
-  //  cout << "TUHKMgen::SetAllParameters() IN" << endl;
-
   fInitialState->fParams.fSqrtS = fHydjetParams.fSqrtS;
   fInitialState->fParams.fAw = fHydjetParams.fAw;
   fInitialState->fParams.fIfb = fHydjetParams.fIfb ;
@@ -415,8 +372,8 @@ void TUHKMgen::SetAllParameters() {
   fInitialState->fParams.fDelta = fHydjetParams.fDelta;
   fInitialState->fParams.fEpsilon = fHydjetParams.fEpsilon;
 
- // fInitialState->fParams.fTime = fHydjetParams.fTime;
   fInitialState->fParams.fWeakDecay = fHydjetParams.fWeakDecay;
+  fInitialState->fParams.fDecay = 1;                   // run decays
   fInitialState->fParams.fCorrS = fHydjetParams.fCorrS;
   fInitialState->fParams.fEtaType = fHydjetParams.fEtaType;
  
index 4f84384..3f1c582 100755 (executable)
@@ -10,7 +10,7 @@
 #include "TGenerator.h"
 #endif
 
-#ifndef INITIALSTATEHYDJET_INCLUDED
+#ifndef INITIALSTATEHYDJET_H
 #include "UHKM/InitialStateHydjet.h"
 #endif
 
@@ -27,7 +27,6 @@ class TUHKMgen : public TGenerator {
  protected:
   InitialStateHydjet *fInitialState;     // HYDJET++ main class which handles the entire Monte-Carlo simulation
   ParticleAllocator fAllocator;       // object which allocates/deallocates memory for the lists of particles
-  List_t fSourceList;             // list holding the initial particles generated by both soft and hard model components  
   List_t fSecondariesList; // list holding the initial particles and the final state particles generated in resonance decays
   Int_t  fNPprim;          // number of primary particles
   Int_t  fNPsec;           // secondary particles
index 9d187c4..4de0e63 100644 (file)
@@ -1,3 +1,9 @@
+// DatabasePDG stores and handles PDG information
+// The PDG particle definitions and decay channels are read
+// in the begining from ASCII files
+// PDG definitions loaded can be selected according to their
+// mass and decay width
+
 /*
   Copyright   : The FASTMC and SPHMC Collaboration
   Author      : Ionut Cristian Arsene 
@@ -11,7 +17,7 @@
 */
 
 
-#ifndef DATABASE_PDG
+#ifndef DATABASEPDG_H
 #include "DatabasePDG.h"
 #endif
 
@@ -33,6 +39,7 @@ DatabasePDG::DatabasePDG():
   fMinimumMass(0.0001),
   fMaximumMass(10.)
 {
+  // Default constructor, initialize members, set input files  
   
   strcpy(fParticleFilename, "particles.data");
   strcpy(fDecayFilename, "tabledecay.txt");
@@ -61,6 +68,8 @@ Bool_t DatabasePDG::LoadData() {
 }
 
 Bool_t DatabasePDG::LoadParticles() {
+  // Read particle definitions from the ascii file  
+
   ifstream particleFile;
   particleFile.open(fParticleFilename);
   if(!particleFile) {
@@ -131,39 +140,41 @@ Bool_t DatabasePDG::LoadParticles() {
 }
 
 Bool_t DatabasePDG::LoadDecays() {
+  // Read the decay channel definitions from the ascii file
+
   ifstream decayFile;
   decayFile.open(fDecayFilename);
   if(!decayFile) {
     cout << "ERROR in DatabasePDG::LoadDecays() : The ASCII file containing the decays list (\""
-         << fDecayFilename << "\") was not found !! Aborting..." << endl;
+         << fDecayFilename << "\") was not found !! Exiting..." << endl;
     return kFALSE;
   }
   
-  Int_t mother_pdg, daughter_pdg[3];
+  Int_t motherPdg, daughterPdg[3];
   Double_t branching;
   
   decayFile.exceptions(ios::failbit);
   while(!decayFile.eof()) {
-    mother_pdg = 0;
-    for(Int_t i=0; i<3; i++) daughter_pdg[i] = 0;
+    motherPdg = 0;
+    for(Int_t i=0; i<3; i++) daughterPdg[i] = 0;
     branching = -1.0;
     try {
-      decayFile >> mother_pdg;
+      decayFile >> motherPdg;
       for(Int_t i=0; i<3; i++) 
-        decayFile >> daughter_pdg[i];
+        decayFile >> daughterPdg[i];
       decayFile >> branching;
     }
     catch (ios::failure const &problem) {
       cout << problem.what() << endl;
       break;
     }
-    if((mother_pdg!=0) && (daughter_pdg[0]!=0) && (branching>=0)) {
+    if((motherPdg!=0) && (daughterPdg[0]!=0) && (branching>=0)) {
       Int_t nDaughters = 0;
       for(Int_t i=0; i<3; i++)
-        if(daughter_pdg[i]!=0)
+        if(daughterPdg[i]!=0)
           nDaughters++;
-      ParticlePDG* particle = GetPDGParticle(mother_pdg);
-      DecayChannel decay(mother_pdg, branching, nDaughters, daughter_pdg);
+      ParticlePDG* particle = GetPDGParticle(motherPdg);
+      DecayChannel decay(motherPdg, branching, nDaughters, daughterPdg);
       particle->AddChannel(&decay);
     }
   }
@@ -176,7 +187,9 @@ Bool_t DatabasePDG::LoadDecays() {
   return kTRUE;
 }
 
-ParticlePDG* DatabasePDG::GetPDGParticleByIndex(Int_t index) {
+ParticlePDG* DatabasePDG::GetPDGParticleByIndex(Int_t index) const {
+  // Return a PDG particle definition based on its index in the particle list 
+
   if(index<0 || index>fNParticles) {
     cout << "Warning in DatabasePDG::GetPDGParticleByIndex(Int_t): Particle index is negative or too big !!" << endl
          << " It must be inside this range: [0, " << fNParticles-1 << "]" << endl
@@ -186,7 +199,11 @@ ParticlePDG* DatabasePDG::GetPDGParticleByIndex(Int_t index) {
   return fParticles[index];
 }
 
-Bool_t DatabasePDG::GetPDGParticleStatusByIndex(Int_t index) {
+Bool_t DatabasePDG::GetPDGParticleStatusByIndex(Int_t index) const {
+  // Return the status of a PDG particle definition based on its index in the particle list
+  // The status is kTRUE when a particle passed the mass, width and charm criteria
+  // and kFALSE otherwise
+
   if(index<0 || index>fNParticles) {
     cout << "Warning in DatabasePDG::GetPDGParticleStatusByIndex(Int_t): Particle index is negative or too big !!" << endl
          << " It must be inside this range: [0, " << fNParticles-1 << "]" << endl
@@ -196,7 +213,10 @@ Bool_t DatabasePDG::GetPDGParticleStatusByIndex(Int_t index) {
   return fStatus[index];
 }
 
-ParticlePDG* DatabasePDG::GetPDGParticle(Int_t pdg) {
+ParticlePDG* DatabasePDG::GetPDGParticle(Int_t pdg) const {
+  // Return a PDG particle definition based on the PDG code (PYTHIA convention) 
+  // If more than 1 definition with the asked PDG code is found then a warning is issued 
+
   Int_t nFindings = 0;
   Int_t firstTimeIndex = 0;
   for(Int_t i=0; i<fNParticles; i++) {
@@ -220,7 +240,11 @@ ParticlePDG* DatabasePDG::GetPDGParticle(Int_t pdg) {
   return 0x0;
 }
 
-Bool_t DatabasePDG::GetPDGParticleStatus(Int_t pdg) {
+Bool_t DatabasePDG::GetPDGParticleStatus(Int_t pdg) const {
+  // Return the status of a PDG particle definition based on its PDG code (PYTHIA convention)
+  // The status is kTRUE when a particle passed the mass, width and charm criteria
+  // and kFALSE otherwise
+
   Int_t nFindings = 0;
   Int_t firstTimeIndex = 0;
   for(Int_t i=0; i<fNParticles; i++) {
@@ -244,7 +268,10 @@ Bool_t DatabasePDG::GetPDGParticleStatus(Int_t pdg) {
   return kFALSE;
 }
 
-ParticlePDG* DatabasePDG::GetPDGParticle(Char_t* name) {
+ParticlePDG* DatabasePDG::GetPDGParticle(Char_t* name) const {
+  // Return a PDG particle definition based on its name
+  // If more than 1 definition with the asked PDG code is found then a warning is issued
+
   Int_t nFindings = 0;
   Int_t firstTimeIndex = 0;
   for(Int_t i=0; i<fNParticles; i++) {
@@ -268,7 +295,11 @@ ParticlePDG* DatabasePDG::GetPDGParticle(Char_t* name) {
   return 0x0;
 }
 
-Bool_t DatabasePDG::GetPDGParticleStatus(Char_t* name) {
+Bool_t DatabasePDG::GetPDGParticleStatus(Char_t* name) const {
+  // Return the status of a PDG particle definition based on its name
+  // The status is kTRUE when a particle passed the mass, width and charm criteria
+  // and kFALSE otherwise
+
   Int_t nFindings = 0;
   Int_t firstTimeIndex = 0;
   for(Int_t i=0; i<fNParticles; i++) {
@@ -292,7 +323,9 @@ Bool_t DatabasePDG::GetPDGParticleStatus(Char_t* name) {
   return kFALSE;
 }
 
-void DatabasePDG::DumpData(Bool_t dumpAll) {
+void DatabasePDG::DumpData(Bool_t dumpAll) const {
+  // Printout all the information in the PDG database
+
   cout << "***********************************************************************************************************" << endl;
   cout << "Dumping all the information contained in the database..." << endl;
   Int_t nDecays = 0;
@@ -346,8 +379,9 @@ void DatabasePDG::DumpData(Bool_t dumpAll) {
   }
 }
 
-Int_t DatabasePDG::CheckImpossibleDecays(Bool_t dump) {
+Int_t DatabasePDG::CheckImpossibleDecays(Bool_t dump) const {
   // Check the database for impossible decays
+
   Int_t nImpossibleDecays = 0;
   for(Int_t currPart=0; currPart<fNParticles; currPart++) {
     if(!fStatus[currPart]) continue;
@@ -394,6 +428,8 @@ Int_t DatabasePDG::CheckImpossibleDecays(Bool_t dump) {
 }
 
 void DatabasePDG::SetUseCharmParticles(Bool_t flag) {
+  // Switch on/off the use of charmed particles
+
   if(fNParticles>0) {
     fUseCharmParticles = flag;
     for(Int_t i=0; i<fNParticles; i++) {
@@ -409,6 +445,7 @@ void DatabasePDG::SetUseCharmParticles(Bool_t flag) {
 }
 
 void DatabasePDG::SetMinimumWidth(Double_t value) {
+  // Set the minimum decay width for the particle definitions to be generated in the soft fireball
   if(fNParticles>0) {
     fMinimumWidth = value;
     for(Int_t i=0; i<fNParticles; i++) {
@@ -424,6 +461,8 @@ void DatabasePDG::SetMinimumWidth(Double_t value) {
 }
 
 void DatabasePDG::SetMaximumWidth(Double_t value) {
+  // Set the maximum decay width for the particle definitions to be generated in the soft fireball
+
   if(fNParticles>0) {
     fMaximumWidth = value;
     for(Int_t i=0; i<fNParticles; i++) {
@@ -439,6 +478,8 @@ void DatabasePDG::SetMaximumWidth(Double_t value) {
 }
 
 void DatabasePDG::SetWidthRange(Double_t min, Double_t max) {
+  // Set the decay width range for the particle definitions to be generated in the soft fireball
+
   if(fNParticles>0) {
     fMinimumWidth = min;
     fMaximumWidth = max;
@@ -457,6 +498,8 @@ void DatabasePDG::SetWidthRange(Double_t min, Double_t max) {
 }
 
 void DatabasePDG::SetMinimumMass(Double_t value) {
+  // Set the minimum mass for the particle definitions to be generated in the soft fireball
+
   if(fNParticles>0) {
     fMinimumMass = value;
     for(Int_t i=0; i<fNParticles; i++) {
@@ -472,6 +515,8 @@ void DatabasePDG::SetMinimumMass(Double_t value) {
 }
 
 void DatabasePDG::SetMaximumMass(Double_t value) {
+  // Set the maximum mass for the particle definitions to be generated in the soft fireball
+
   if(fNParticles>0) {
     fMaximumMass = value;
     for(Int_t i=0; i<fNParticles; i++) {
@@ -487,6 +532,8 @@ void DatabasePDG::SetMaximumMass(Double_t value) {
 }
 
 void DatabasePDG::SetMassRange(Double_t min, Double_t max) {
+  // Set the mass range for the particle definitions to be generated in the soft fireball
+
   if(fNParticles>0) {
     fMinimumMass = min;
     fMaximumMass = max;
@@ -505,6 +552,8 @@ void DatabasePDG::SetMassRange(Double_t min, Double_t max) {
 }
 
 void DatabasePDG::SortParticles() {
+  // Sort the particle list so that those with kTRUE status will be always on top of the list
+
   if(fNParticles<2) {
     cout << "Warning in DatabasePDG::SortParticles() : No particles to sort. Load data first!!" << endl;
     return;
@@ -536,7 +585,11 @@ void DatabasePDG::SortParticles() {
   return;
 }
 
-Int_t DatabasePDG::GetNParticles(Bool_t all) {
+Int_t DatabasePDG::GetNParticles(Bool_t all) const {
+  // Return the number of particle definitions in the database
+  // If all is kTRUE then return number of all particle
+  // If all is kFALSE then return the number of good status (kTRUE) particles
+
   if(all)
     return fNParticles;
 
@@ -547,6 +600,9 @@ Int_t DatabasePDG::GetNParticles(Bool_t all) {
 }
 
 void DatabasePDG::UseThisListOfParticles(Char_t *filename, Bool_t exclusive) {
+  // Read a list of PDG codes from the file "filename" and mark them with good status (kTRUE)
+  // while all the other will be marked kFALSE (only if exclusive = kTRUE)     
+
   if(fNParticles<1) {
     cout << "Error in DatabasePDG::UseThisListOfParticles(Char_t*, Bool_t) : You must load the data before calling this function!!" << endl;
     return;
@@ -556,7 +612,7 @@ void DatabasePDG::UseThisListOfParticles(Char_t *filename, Bool_t exclusive) {
   listFile.open(filename);
   if(!listFile) {
     cout << "ERROR in DatabasePDG::UseThisListOfParticles(Char_t*, Bool_t) : The ASCII file containing the PDG codes list (\""
-         << filename << "\") was not found !! Aborting..." << endl;
+         << filename << "\") was not found !! Exiting..." << endl;
     return;
   }
 
@@ -603,7 +659,9 @@ void DatabasePDG::UseThisListOfParticles(Char_t *filename, Bool_t exclusive) {
   return;
 }
 
-Bool_t DatabasePDG::IsChannelAllowed(DecayChannel *channel, Double_t motherMass) {
+Bool_t DatabasePDG::IsChannelAllowed(DecayChannel *channel, Double_t motherMass) const {
+  // Check if the decay channel "channel" is allowed by using the mother particle mass "motherMass"
+
   Double_t daughtersSumMass = 0.0;
   for(Int_t i=0; i<channel->GetNDaughters(); i++)
     daughtersSumMass += GetPDGParticle(channel->GetDaughterPDG(i))->GetMass();
@@ -612,7 +670,9 @@ Bool_t DatabasePDG::IsChannelAllowed(DecayChannel *channel, Double_t motherMass)
   return kFALSE;
 }
 
-Int_t DatabasePDG::GetNAllowedChannels(ParticlePDG *particle, Double_t motherMass) {
+Int_t DatabasePDG::GetNAllowedChannels(ParticlePDG *particle, Double_t motherMass) const {
+  // Check how many decay channels are allowed for a given particle definition at a given mass
+
   Int_t nAllowedChannels = 0;
   for(Int_t i=0; i<particle->GetNDecayChannels(); i++)
     nAllowedChannels += (IsChannelAllowed(particle->GetDecayChannel(i), motherMass) ? 1:0);
index f0bdda0..89efc3d 100644 (file)
@@ -1,3 +1,7 @@
+// DatabasePDG stores and handles PDG information
+// The PDG particle definitions and decay channels are read
+// in the begining from ASCII files
+
 /*
   Copyright   : The FASTMC and SPHMC Collaboration
   Author      : Ionut Cristian Arsene 
   SHARE (Computer Physics Communications 167 229 (2005)) collaborations.
 */
 
-#ifndef DATABASE_PDG
-#define DATABASE_PDG
+#ifndef DATABASEPDG_H
+#define DATABASEPDG_H
 
 #include "Rtypes.h"
 #ifndef PARTICLE_PDG
 #include "ParticlePDG.h"
 #endif
 
-const Int_t kMaxParticles = 500;
+const Int_t kMaxParticles = 1000;
 
 class DatabasePDG {
  public:
@@ -54,24 +58,24 @@ class DatabasePDG {
   
   Char_t* GetParticleFilename() {return fParticleFilename;}
   Char_t* GetDecayFilename() {return fDecayFilename;}
-  Int_t GetNParticles(Bool_t all = kFALSE);      // true - no. of all particles; false - no. of good status particles
-  ParticlePDG* GetPDGParticleByIndex(Int_t index);
-  Bool_t GetPDGParticleStatusByIndex(Int_t index);
-  ParticlePDG* GetPDGParticle(Int_t pdg);
-  Bool_t GetPDGParticleStatus(Int_t pdg);
-  ParticlePDG* GetPDGParticle(Char_t *name);
-  Bool_t GetPDGParticleStatus(Char_t *name);
-  Bool_t GetUseCharmParticles() {return fUseCharmParticles;};
-  Double_t GetMinimumWidth() {return fMinimumWidth;};
-  Double_t GetMaximumWidth() {return fMaximumWidth;};
-  Double_t GetMinimumMass() {return fMinimumMass;};
-  Double_t GetMaximumMass() {return fMaximumMass;};
-  void DumpData(Bool_t dumpAll = kFALSE); // print the PDG information in the console
-  Int_t CheckImpossibleDecays(Bool_t dump = kFALSE);   // print all impossible decays included in the database
-  Bool_t IsChannelAllowed(DecayChannel *channel, Double_t motherMass);
-  Int_t GetNAllowedChannels(ParticlePDG *particle, Double_t motherMass);
+  Int_t GetNParticles(Bool_t all = kFALSE) const;      // true - no. of all particles; false - no. of good status particles
+  ParticlePDG* GetPDGParticleByIndex(Int_t index) const;
+  Bool_t GetPDGParticleStatusByIndex(Int_t index) const;
+  ParticlePDG* GetPDGParticle(Int_t pdg) const;
+  Bool_t GetPDGParticleStatus(Int_t pdg) const;
+  ParticlePDG* GetPDGParticle(Char_t *name) const;
+  Bool_t GetPDGParticleStatus(Char_t *name) const; 
+  Bool_t GetUseCharmParticles() const {return fUseCharmParticles;};
+  Double_t GetMinimumWidth() const {return fMinimumWidth;};
+  Double_t GetMaximumWidth() const {return fMaximumWidth;};
+  Double_t GetMinimumMass() const {return fMinimumMass;};
+  Double_t GetMaximumMass() const {return fMaximumMass;};
+  void DumpData(Bool_t dumpAll = kFALSE) const; // print the PDG information in the console
+  Int_t CheckImpossibleDecays(Bool_t dump = kFALSE) const;   // print all impossible decays included in the database
+  Bool_t IsChannelAllowed(DecayChannel *channel, Double_t motherMass) const;
+  Int_t GetNAllowedChannels(ParticlePDG *particle, Double_t motherMass) const;
   void SetStable(Int_t pdg, Bool_t value) {GetPDGParticle(pdg)->SetStable(value);}
-  Bool_t GetStableStatus(Int_t pdg) {return GetPDGParticle(pdg)->GetStableStatus();}
+  Bool_t GetStableStatus(Int_t pdg) const {return GetPDGParticle(pdg)->GetStableStatus();}
 
  private:
   DatabasePDG(const DatabasePDG&);
index 4784e64..17fb3d0 100644 (file)
@@ -12,7 +12,7 @@
 #ifndef PARTICLE_PDG
 #include "ParticlePDG.h"
 #endif
-#ifndef DATABASE_PDG
+#ifndef DATABASEPDG_H
 #include "DatabasePDG.h"
 #endif
 
index cc2417d..2e04bfd 100644 (file)
@@ -17,7 +17,7 @@
 #include <TError.h>
 #include <TMath.h>
 
-#ifndef DATABASE_PDG
+#ifndef DATABASEPDG_H
 #include "DatabasePDG.h"
 #endif
 #ifndef PARTICLE_PDG
@@ -64,14 +64,6 @@ extern "C" void mydelta_();
 extern SERVICEEVCommon SERVICEEV;
 
 void Decay(List_t &output, Particle &parent, ParticleAllocator &allocator, DatabasePDG* database) {
-  // check if the parent particle has been decayed already
-  //  std::cout << "HadronDecayer::Decay() IN" << std::endl;
-  Int_t daughters = parent.GetNDaughters();
-
-  //cout<<"in Decay pdg"<< parent.Def()->GetPDG()<<"daughters "<<daughters<<endl;
-
-  if(daughters>0)  // particle decayed already
-    return;
 
   // Get the PDG properties of the particle
   ParticlePDG *pDef = parent.Def();
@@ -88,7 +80,7 @@ void Decay(List_t &output, Particle &parent, ParticleAllocator &allocator, Datab
 
   // get the PDG mass of the specie  
   Double_t PDGmass = pDef->GetMass();
-  int ComprCodePyth=0;
+  Int_t ComprCodePyth=0;
   Float_t Delta =0;
 
   Bool_t success = kFALSE;
@@ -102,9 +94,9 @@ void Decay(List_t &output, Particle &parent, ParticleAllocator &allocator, Datab
       std::cout << "               Will be left undecayed ... check it out!" << std::endl;
       return;
     }
+    
     // get a random mass using the Breit Wigner distribution 
-    Double_t BWmass = gRandom->BreitWigner(PDGmass, pDef->GetWidth()); 
-    //!!!!    
+    Double_t BWmass = gRandom->BreitWigner(PDGmass, pDef->GetWidth());     
     //      BWmass = PDGmass;
     // Try to cut the Breit Wigner tail of the particle using the cuts from pythia
     // The Delta variable is obtained from pythia based on the specie
@@ -132,26 +124,18 @@ void Decay(List_t &output, Particle &parent, ParticleAllocator &allocator, Datab
       Delta=0.0;
     }
 
-    //    std::cout << "HadronDecayer::Decay() BWmass = " << BWmass << std::endl;
-    //    std::cout << "HadronDecayer::Decay() Delta = " << Delta << std::endl;
     //for particles from PYTHIA table only, if the BW mass is outside the cut range then quit this iteration and generate another BW mass
     if(ComprCodePyth!=0 && Delta>0 && (BWmass<PDGmass-Delta || BWmass>PDGmass+Delta)){
-      //      std::cout<<"encoding"<<encoding<<"delta"<<Delta<<"width "<<pDef->GetWidth()<<"mass"<<BWmass<<std::endl;
       iterations++;
-      //      std::cout << "HadronDecayer::Decay() BWmass outside cut, try again" << std::endl;
       continue;
     }    
-    //----    
-    
-    //    if(BWmass>5)
-    //      std::cout<<" > 5 encoding"<<encoding<<" pdgmass "<<PDGmass<<" delta "<<Delta<<"width "<<pDef->GetWidth()<<" mass "<<BWmass<<"CC"<<ComprCodePyth<<std::endl;
     
     // check how many decay channels are allowed with the generated mass
     Int_t nAllowedChannels = database->GetNAllowedChannels(pDef, BWmass);
     // if no decay channels are posible with this mass, then generate another BW mass
     if(nAllowedChannels==0) {    
       iterations++;
-      //      std::cout << "HadronDecayer::Decay() no decays allowed at this BW mass" << std::endl;
+            std::cout << "HadronDecayer::Decay() no decays allowed at this BW mass" << std::endl;
       continue;
     }
 
@@ -166,14 +150,10 @@ void Decay(List_t &output, Particle &parent, ParticleAllocator &allocator, Datab
     Int_t chosenChannel = 1000;
     Bool_t found = kFALSE;
     Int_t channelIterations = 0;
-    //    std::cout << "HadronDecayer::Decay() randValue = " << randValue  << std::endl;
     while(!found) {
-      //      std::cout << "HadronDecayer::Decay() channel iteration #" << channelIterations << std::endl;
-      for(Int_t nChannel = 0; nChannel < nDecayChannel; ++nChannel) {
+      for(Int_t nChannel = 0; nChannel < nDecayChannel; nChannel++) {
        randValue -= pDef->GetDecayChannel(nChannel)->GetBranching();
-       //      std::cout << "HadronDecayer::Decay() channel #" << nChannel << " randValue = " << randValue << std::endl;
        if(randValue <= 0. && database->IsChannelAllowed(pDef->GetDecayChannel(nChannel), BWmass)) {
-         //      std::cout << "HadronDecayer::Decay() channel found" << std::endl;
          chosenChannel = nChannel;
          found = kTRUE;
          break;
@@ -213,14 +193,12 @@ void Decay(List_t &output, Particle &parent, ParticleAllocator &allocator, Datab
       p1.SetLastMotherDecayMom(parentBW.Mom());
       p1.SetType(NB);
 
-      // link the parent and daughters trough their indexes in the list
-      Int_t parentIndex = -1;
-      if(parent.GetMother()==-1) parentIndex = parent.SetIndex();   // parents which are primaries don't have yet an index
-      else parentIndex = parent.GetIndex();                         // parents which are secondaries have an index
+      // add the daughter to the list of secondaries
+      Int_t parentIndex = parent.GetIndex();
       Int_t p1Index = p1.SetIndex();                           // set the daughter index
       p1.SetMother(parentIndex);                               // set the mother index for this daughter 
-      parent.SetDaughter(p1Index);                             // set p1 as daughter to the parent 
-      if(parent.GetMother()==-1) allocator.AddParticle(parent, output);   // add it only if its a primary particle
+      parent.SetFirstDaughterIndex(p1Index);
+      parent.SetLastDaughterIndex(p1Index);
       allocator.AddParticle(p1, output);
       success = kTRUE;  
     }
@@ -246,7 +224,6 @@ void Decay(List_t &output, Particle &parent, ParticleAllocator &allocator, Datab
       p2.SetLastMotherPdg(parentBW.Encoding());
       p2.SetLastMotherDecayCoor(parentBW.Pos());
       p2.SetLastMotherDecayMom(parentBW.Mom());
-      // std::cout<<"2d NB="<<NB<<std::endl;
       //set to daughters the same type as has mother
       p1.SetType(NB);
       p2.SetType(NB);
@@ -259,35 +236,22 @@ void Decay(List_t &output, Particle &parent, ParticleAllocator &allocator, Datab
                                    (parentBW.Mom().E()-p1.Mom().E()-p2.Mom().E())*(parentBW.Mom().E()-p1.Mom().E()-p2.Mom().E()));
       // if deltaS is too big then repeat the kinematic procedure
  
       if(deltaS>0.001) {
-
-       //      cout << "2-body decay kinematic check in lab system: " << pDef->GetPDG() << " --> " << p1.Encoding() << " + " << p2.Encoding() << endl;
-       //      cout << "Mother    (e,px,py,pz): " << parentBW.Mom().E() << "\t" << parentBW.Mom().X() << "\t" << parentBW.Mom().Y() << "\t" << parentBW.Mom().Z() << endl;
-       //      cout << "Mother    (x,y,z,t): " << parentBW.Pos().X() << "\t" << parentBW.Pos().Y() << "\t" << parentBW.Pos().Z() << "\t" << parentBW.Pos().T() << endl;
-
-       //      cout << "Daughter1 (e,px,py,pz): " << p1.Mom().E() << "\t" << p1.Mom().X() << "\t" << p1.Mom().Y() << "\t" << p1.Mom().Z() << endl;
-       //      cout << "Daughter2 (e,px,py,pz): " << p2.Mom().E() << "\t" << p2.Mom().X() << "\t" << p2.Mom().Y() << "\t" << p2.Mom().Z() << endl;     
-       //      cout << "2-body decay delta(sqrtS) = " << deltaS << endl;
-       //      cout << "Repeating the decay algorithm ..." << endl;
-
        iterations++;
        continue;
       }
       // push particles to the list of secondaries
-      Int_t parentIndex = -1;
-      if(parent.GetMother()==-1) parentIndex = parent.SetIndex();   // parents which are primaries don't have yet an index
-      else parentIndex = parent.GetIndex();                         // parents which are secondaries have an index
+      Int_t parentIndex = parent.GetIndex();
       p1.SetIndex(); 
       p2.SetIndex();
       p1.SetMother(parentIndex); 
       p2.SetMother(parentIndex);
-      parent.SetDaughter(p1.GetIndex());
-      parent.SetDaughter(p2.GetIndex());
-      if(parent.GetMother()==-1) allocator.AddParticle(parent, output);      // add it only if its a primary
+      parent.SetFirstDaughterIndex(p1.GetIndex());
+      parent.SetLastDaughterIndex(p2.GetIndex());
       allocator.AddParticle(p1, output);
       allocator.AddParticle(p2, output);
       success = kTRUE;
+    
     }
 
     // third case: three daughter particle
@@ -381,10 +345,7 @@ void Decay(List_t &output, Particle &parent, ParticleAllocator &allocator, Datab
       p1.SetType(NB);
       p2.SetType(NB);
       p3.SetType(NB);
- //      std::cout<<"3d NB="<<NB<<std::endl;
-
 
-            
       // energy conservation check in the lab system
       Double_t deltaS = TMath::Sqrt((parentBW.Mom().X()-p1.Mom().X()-p2.Mom().X()-p3.Mom().X())*(parentBW.Mom().X()-p1.Mom().X()-p2.Mom().X()-p3.Mom().X()) +
                                    (parentBW.Mom().Y()-p1.Mom().Y()-p2.Mom().Y()-p3.Mom().Y())*(parentBW.Mom().Y()-p1.Mom().Y()-p2.Mom().Y()-p3.Mom().Y()) +
@@ -392,7 +353,6 @@ void Decay(List_t &output, Particle &parent, ParticleAllocator &allocator, Datab
                                    (parentBW.Mom().E()-p1.Mom().E()-p2.Mom().E()-p3.Mom().E())*(parentBW.Mom().E()-p1.Mom().E()-p2.Mom().E()-p3.Mom().E()));
       // if deltaS is too big then repeat the kinematic procedure
       if(deltaS>0.001) {
-
        //      cout << "3-body decay kinematic check in lab system: " << pDef->GetPDG() << " --> " << p1.Encoding() << " + " << p2.Encoding() << " + " << p3.Encoding() << endl;
        //      cout << "Mother    (e,px,py,pz): " << parentBW.Mom().E() << "\t" << parentBW.Mom().X() << "\t" << parentBW.Mom().Y() << "\t" << parentBW.Mom().Z() << endl;
        //      cout << "Daughter1 (e,px,py,pz): " << p1.Mom().E() << "\t" << p1.Mom().X() << "\t" << p1.Mom().Y() << "\t" << p1.Mom().Z() << endl;
@@ -405,19 +365,15 @@ void Decay(List_t &output, Particle &parent, ParticleAllocator &allocator, Datab
        continue;
       }
 
-      Int_t parentIndex = -1;
-      if(parent.GetMother()==-1) parentIndex = parent.SetIndex();   // parents which are primaries don't have yet an index
-      else parentIndex = parent.GetIndex();                         // parents which are secondaries have an index
+      Int_t parentIndex = parent.GetIndex();
       p1.SetIndex();
       p2.SetIndex();
       p3.SetIndex();
       p1.SetMother(parentIndex); 
       p2.SetMother(parentIndex);
       p3.SetMother(parentIndex);
-      parent.SetDaughter(p1.GetIndex());
-      parent.SetDaughter(p2.GetIndex());
-      parent.SetDaughter(p3.GetIndex());
-      if(parent.GetMother()==-1) allocator.AddParticle(parent, output);    // add it only if its a primary
+      parent.SetFirstDaughterIndex(p1.GetIndex());
+      parent.SetLastDaughterIndex(p3.GetIndex());
       allocator.AddParticle(p1, output);
       allocator.AddParticle(p2, output);
       allocator.AddParticle(p3, output);
index 950cd12..0630039 100644 (file)
@@ -8,7 +8,7 @@
                            November. 2, 2005                                
 
 */
-#ifndef DATABASE_PDG
+#ifndef DATABASEPDG_H
 #include "DatabasePDG.h"
 #endif
 #ifndef PARTICLE_INCLUDED
index b5a08cc..0e55891 100644 (file)
 #include <iostream> 
 #include <fstream>
 using namespace std;
+
 // aic(2008/08/08): define the fLastIndex static variable 
-Int_t Particle::fLastIndex;
+//Int_t Particle::fLastIndex;
+
+void InitialState::Evolve(List_t &secondaries, ParticleAllocator &allocator, Double_t weakDecayLimit) {
+  // Particle indexes are set for primaries already from InitialStateHydjet::Initialize()
 
-void InitialState::Evolve(List_t &source, List_t &secondaries, ParticleAllocator &allocator, Double_t weakDecayLimit) {
-  cout << "InitialState::Evolve()  IN" << endl;
+  // particle list iterators
   LPIT_t it;
   LPIT_t e;
 
-  // Initialize the fLastIndex to -1
-  it = source.begin();
-  it->InitIndexing();
+  // Particle indexes are set for primaries already from InitialStateHydjet::Initialize()
 
-  //  cout << "InitialState::Evolve()  Loop over primaries" << endl;
-  // Loop over source list (primary particles) 
-  for(it = source.begin(), e = source.end(); it != e; ++it) {
-    //    cout << "InitialState::Evolve()  particle(in primary loop) code = " << it->Encoding() << endl;
-    //    if(GetPDGParticleStableStatus(it->Encoding())) {
-      //      cout << "InitialState::Evolve() its a stable marked particle, decay time should be 0" << endl;
-      //    }
-    // calculate the decay time of the particle
-    Double_t decayTime = GetDecayTime(*it, weakDecayLimit);
-    //    cout << "InitialState::Evolve()  decay time = " << decayTime << endl;
-    TVector3 shift(it->Mom().BoostVector());
-    shift *= decayTime;
-    it->SetLastInterTime(it->T() + decayTime);
+  // Decay loop
+  // Note that the decay products are always added at the end of list so the unstable products are
+  // decayed when the iterator reaches them (used for cascade decays)
 
-    // if the particle is unstable run the decayer 
-    if(decayTime > 0.) {
-      //      cout << "InitialState::Evolve()  perform decay..." << endl;
-      it->Pos(it->Pos() += TLorentzVector(shift, 0.));
-      it->T(it->T() + decayTime);
-      // perform the decay procedure
-      // the primaries are added to the list of secondaries inside Decay()
-      Decay(secondaries, *it, allocator, fDatabase);  
-    }
-    // else just add the primary to the list of particles
-    else {
-      //      cout << "InitialState::Evolve()  stable particle (just add it to the list of final particles)" << endl;
-      it->SetIndex();
-      allocator.AddParticle(*it, secondaries);
-    }
-  }
+  it = secondaries.begin();
 
-  //  cout << "InitialState::Evolve()  Loop over secondaries" << endl;
-  // Loop over the secondaries list and decay the cascades
-  for(it = secondaries.begin(), e = secondaries.end(); it != e;) {
-    //    cout << "InitialState::Evolve()  particle(in secondaries loop) code " << it->Encoding() << endl;
-    if(GetPDGParticleStableStatus(it->Encoding())) {
-      //      cout << "InitialState::Evolve()  its a stable particle " << it->Encoding()
-      //          << "(already added to final list)"  << endl;
-      ++it;
-      continue;
-    }
-    // the primaries are already decayed, so just ignore it
-    if(it->GetMother()==-1) {  
-      //      cout << "InitialState::Evolve()  its a primary (already decayed)" << endl;
-      ++it;
+  for(it = secondaries.begin(), e = secondaries.end(); it != e; it++) {
+    // if the decay procedure was applied already skip ... (e.g. particles from pythia history information)
+    if(it->GetDecayed()) {
       continue;
     }
+
+    // generate the decay time; if particle is stable or set to stable decayTime=0
     Double_t decayTime = GetDecayTime(*it, weakDecayLimit);
-    //    cout << "InitialState::Evolve()  decay time = " << decayTime << endl;
     it->SetLastInterTime(it->T() + decayTime);
     TVector3 shift(it->Mom().BoostVector());
     shift *= decayTime; 
-    
-    // if the particle is unstable run the decays
+    it->SetDecayed();
+
+    // if decayTime>0 then apply the decay procedure (only 2 or 3 body decays)
     if(decayTime > 0.) {
-      //      cout << "InitialState::Evolve()  perform decay..." << endl;
       it->Pos(it->Pos() += TLorentzVector(shift, 0.));
       it->T(it->T() + decayTime);
-      // The decay products are added to the list inside the Decay() function
-      // The mother particle is not removed from list but can be distinguished 
-      // since the daughter indexes will be initialized
       Decay(secondaries, *it, allocator, fDatabase);
-      ++it;
-    }
-    // else just continue
-    else {
-      //      cout << "InitialState::Evolve() no decay, continue to next particle" << endl;
-      ++it;
     }
+    // if particle is stable just continue
+  }
+
+  it = secondaries.begin();
+  Int_t npart = it->GetLastIndex();
+  cout << "particles generated = " << npart << endl;
+  for(Int_t count=0; count<npart; count++, it++) {
+//    cout << "=====================================================" << endl;
+//    cout << "InitialState::Evolve() particle count = " << count << endl;
+//    cout << "InitialState::Evolve() particle pdg = " << it->Encoding() << endl;
+//    cout << "InitialState::Evolve() index = " << it->GetIndex() << endl;
+//    cout << "InitialState::Evolve() mother index = " << it->GetMother() << endl;
+//    cout << "InitialState::Evolve() mother pdg = " << it->GetLastMotherPdg() << endl;
+//    cout << "InitialState::Evolve() n daughters = " << it->GetNDaughters() << endl;
+//    cout << "InitialState::Evolve() d (first, last) = " << it->GetFirstDaughterIndex() << ", " 
+//      << it->GetLastDaughterIndex() << endl;
   }
-  //  cout << "InitialState::Evolve()  OUT"<< endl;
 }
index 4110fc6..f52db7b 100644 (file)
@@ -12,7 +12,7 @@
 #ifndef PARTICLE_INCLUDED
 #include "Particle.h"
 #endif
-#ifndef DATABASE_PDG
+#ifndef DATABASEPDG_H
 #include "DatabasePDG.h"
 #endif
 
@@ -47,11 +47,12 @@ class InitialState {
   virtual void Initialize(List_t &source, ParticleAllocator &allocator) = 0;
   virtual Bool_t ReadParams() = 0;
   virtual Bool_t MultIni() = 0;
-  virtual Double_t GetTime() = 0;
+  virtual Bool_t RunDecays() = 0;
   virtual Int_t GetNev() = 0;
   virtual Double_t GetWeakDecayLimit() = 0;
     
-  virtual void Evolve(List_t &source, List_t &secondaries, ParticleAllocator &allocator, Double_t weakDecayLimit);
+  //  virtual void Evolve(List_t &source, List_t &secondaries, ParticleAllocator &allocator, Double_t weakDecayLimit);
+  virtual void Evolve(List_t &secondaries, ParticleAllocator &allocator, Double_t weakDecayLimit);
  protected:
    DatabasePDG *fDatabase;
  private:
index a16f3fa..c2c50d3 100644 (file)
@@ -1,3 +1,7 @@
+//expanding localy equilibated fireball with volume hadron radiation
+//thermal part: Blast wave model, Bjorken-like parametrization
+//hyght-pt: PYTHIA + jet quenching model PYQUEN
+
 /*                                                                           
          HYDJET++ 
          version 1.0:  
                      
 */
 
-
-//expanding localy equilibated fireball with volume hadron radiation
-//thermal part: Blast wave model, Bjorken-like parametrization
-//hyght-pt: PYTHIA + jet quenching model PYQUEN
-
-
 #include <TLorentzVector.h>
 #include <TVector3.h>
-#include <TRandom.h>
 #include <TMath.h>
 
-#ifndef INITIALSTATEHYDJET_INCLUDED
+#ifndef INITIALSTATEHYDJET_H
 #include "InitialStateHydjet.h"
 #endif
 #ifndef RANDARRAYFUNCTION_INCLUDED
 #include "RandArrayFunction.h"
 #endif
-#ifndef HADRONDECAYER_INCLUDED
-#include "HadronDecayer.h"
-#endif
 #ifndef GRANDCANONICAL_INCLUDED
 #include "GrandCanonical.h"
 #endif
 #ifndef NAStrangePotential_h
 #include "StrangePotential.h"
 #endif
-#ifndef NAEquationSolver_h
-#include "EquationSolver.h"
-#endif
 #ifndef PARTICLE_INCLUDED
 #include "Particle.h"
 #endif
 #ifndef PARTICLE_PDG
 #include "ParticlePDG.h"
 #endif
-#ifndef UKUTILITY_INCLUDED
-#include "UKUtility.h"
-#endif
 #include <iostream> 
 #include <fstream>
 #include "HYJET_COMMONS.h"
 extern "C" void  hyevnt_();
 extern "C" void  myini_();
-//extern "C" void  hyinit_();
 extern HYIPARCommon HYIPAR;
 extern HYFPARCommon HYFPAR;
 extern HYJPARCommon HYJPAR;
@@ -74,35 +61,30 @@ extern SERVICECommon SERVICE;
 using std::cout;
 using std::endl;
 
+class ParticleAllocator;
+class TRandom3;
+
+// declaration of the static member fLastIndex
+Int_t Particle::fLastIndex;
+
 void InitialStateHydjet::Initialize(List_t &source, ParticleAllocator & allocator) {
   // Generate initial particles from the soft and hard components
+
+  // Initialize the static "last index variable"
+  Particle::InitIndexing(); 
+
   //----- high-pt part------------------------------
   TLorentzVector partJMom, partJPos, zeroVec;
-  HYJPAR.ishad  = fParams.fIshad;
-  PYQPAR.T0     = fParams.fT0;
-  PYQPAR.tau0   = fParams.fTau0;
-  PYQPAR.nf     = fParams.fNf;
-  PYQPAR.ienglu = fParams.fIenglu;
-  PYQPAR.ianglu = fParams.fIanglu;
 
   // run a HYDJET event
   hyevnt_(); 
    
-  //get number of particles in jets
-  //Int_t numbJetPart = HYPART.njp;  
-  //Double_t  Bgen = HYFPAR.bgen;
-  //Int_t  Njet = HYJPAR.njet;
-  //Int_t  Nbcol = HYFPAR.nbcol;
-
   if(fParams.fNhsel != 0) {   
+    //get number of particles in jets
     Int_t numbJetPart = HYPART.njp;
-    for(Int_t i = 0; i <numbJetPart; ++i) {
-      Int_t pdg = int(HYPART.ppart[i][1]);
-      //      if(pdg==310) pdg=311; //Kos Kol 130 we have no in the table, we do not put its in the list, the same for e,mu  
-      //      if(pdg==130 || pdg==-130) continue;
-      //      if(pdg==311 || pdg==-311 || pdg=130 || pdg==-130 || pdg==310 || pdg==-310)
-      //       cout << "pdg = " << pdg << endl;
+
+    for(Int_t i = 0; i <numbJetPart; i++) {
+      Int_t pdg = Int_t(HYPART.ppart[i][1]);
       Double_t px = HYPART.ppart[i][2];
       Double_t py = HYPART.ppart[i][3];
       Double_t pz = HYPART.ppart[i][4];
@@ -111,29 +93,22 @@ void InitialStateHydjet::Initialize(List_t &source, ParticleAllocator & allocato
       Double_t vy = HYPART.ppart[i][7];
       Double_t vz = HYPART.ppart[i][8];
       Double_t vt = HYPART.ppart[i][9];    
-      partJMom.SetXYZT(px, py, pz, e); 
-      partJPos.SetXYZT(vx, vy, vz, vt);
-      //std::cout<<" --InitialStateHydjet  pdg "<<pdg<<" px "<<px<<" py "<<py<<" pz "<<pz<<" e"<<e<<std::endl;
-      //        std::cout<<" vx in fm "<<vx<<" vy "<<vy<<" vz "<<vz<<"vt"<<vt<<std::endl;
       ParticlePDG *partDef = fDatabase->GetPDGParticle(pdg);
-      //      cout<<"after ParticlePDG"<<pdg<<" partDef"<<partDef<<endl;
-      
       Int_t type =1;                //from jet
-      if(partDef)
-       allocator.AddParticle(Particle(partDef, partJPos, partJMom, 0, 0,type,-1, zeroVec, zeroVec), source);
-      //      else 
-      //       cout << "particle not added in output (no definition in DatabasePDG)" << endl;
+      if(partDef) {
+       partJMom.SetXYZT(px, py, pz, e);
+       partJPos.SetXYZT(vx, vy, vz, vt);
+       Particle *particle=new Particle(partDef, partJPos, partJMom, 0, 0, type, -1, zeroVec, zeroVec);
+       particle->SetIndex();
+       allocator.AddParticle(*particle, source);
+       delete particle;
+      }
     }
   }       //nhsel !=0 not only hydro!             
 
-  //  std::cout<<"in InitialStateHydjet::Initialize 2"<<std::endl;
-
          
   //----------HYDRO part------------------------------------------------
   if(fParams.fNhsel < 3) {
-   
-    //    cout<<"in HYDRO part 0"<<endl;
-  
     const Double_t  weightMax = 2*TMath::CosH(fParams.fUmax);
     const Int_t nBins = 100;
     Double_t probList[nBins];
@@ -146,8 +121,6 @@ void InitialStateHydjet::Initialize(List_t &source, ParticleAllocator & allocato
     const Double_t eMax = 5.;  
     //-------------------------------------
     // get impact parameter    
-    //    Double_t impactParameter = HYFPAR.bgen;
-
     
     //effective volume for central     
     double dYl= 2 * fParams.fYlmax; //uniform distr. [-Ylmax; Ylmax]  
@@ -181,7 +154,6 @@ void InitialStateHydjet::Initialize(List_t &source, ParticleAllocator & allocato
        }
        //no charm now !
        if(partDef->GetCharmQNumber()!=0 || partDef->GetCharmAQNumber()!=0){
-         //      cout<<"charmed particle, don't use now ! "<<encoding<<endl;
          continue;
        }
 
@@ -279,28 +251,25 @@ void InitialStateHydjet::Initialize(List_t &source, ParticleAllocator & allocato
            weight = (n1 * p0) /e;  // weight for rdr gammar: weight = (n1 * p0) / n1[3] / e; 
          } while(yy >= weight); 
        
-         //  if(abs(z0)>1000 || abs(x0)>1000)std::cout<<"====== etaF==== "<<etaF<<std::endl;
           partMom.SetXYZT(px0, py0, pz0, e);
          partPos.SetXYZT(x0, y0, z0, t0);
          partMom.Boost(vec3);
          Int_t type =0; //hydro
 
-         //           cout<<"before AddParticle in HYDRO part"<<endl;
-         allocator.AddParticle(Particle(partDef, partPos, partMom, 0., 0, type, -1, zeroVec, zeroVec), source);
-         
+         Particle *particle=new Particle(partDef, partPos, partMom, 0., 0, type, -1, zeroVec, zeroVec);
+         particle->SetIndex();
+         allocator.AddParticle(*particle, source);
+         delete particle;
         } //nhsel==4 , no hydro part
       }
     } 
   }
   
-  //    std::cout<<"in InitialStateHydjet::Initialize OUT"<<std::endl;
-
 }
 
 Bool_t InitialStateHydjet::ReadParams() {     
   // Read parameters from an input file in ascii 
  
-  //  std::cout<<"\nWelcome to Hydjet++ hadron freezeout generator!\nInput: \n\n"; 
   Float_t par[200] = {0.};
   Int_t i = 0; 
   std::string s(40,' '); 
@@ -433,7 +402,6 @@ Bool_t InitialStateHydjet::MultIni() {
     //compute strangeness potential
     if(fParams.fMuB > 0.01)
       fParams.fMuS = psp->CalculateStrangePotential();
-    //    cout << "fMuS = " << fParams.fMuS << endl;  
     
     //if user choose fYlmax larger then allowed by kinematics at the specified beam energy sqrt(s)     
     if(fParams.fYlmax > TMath::Log(fParams.fSqrtS/0.94)){
@@ -445,8 +413,6 @@ Bool_t InitialStateHydjet::MultIni() {
     if(fParams.fCorrS <= 0.) {
       //see F. Becattini, J. Mannien, M. Gazdzicki, Phys Rev. C73 044905 (2006)
       fParams.fCorrS = 1. - 0.386* TMath::Exp(-1.23*fParams.fT/fParams.fMuB);
-      //      std::cout<<"The phenomenological f-la F. Becattini et al. PRC73 044905 (2006) for CorrS was used." << std::endl;
-      //      std::cout<<"Strangeness suppression parameter = "<<fParams.fCorrS << std::endl;
       
     }
     std::cout<<"The phenomenological f-la J. Cleymans et al. PRC73 034905 (2006) for Tch mu_B was used." << std::endl;
@@ -461,11 +427,6 @@ Bool_t InitialStateHydjet::MultIni() {
   }
   
   
-  //  std::cout<<"Used eta_max = "<<fParams.fYlmax<<  std::endl;
-  //  std::cout<<"maximal allowed eta_max TMath::Log(fParams.fSqrtS/0.94)=  "<<TMath::Log(fParams.fSqrtS/0.94)<<std::endl;
-  
-  
-  
   //initialisation of high-pt part 
   
   HYJPAR.nhsel = fParams.fNhsel;
@@ -492,33 +453,22 @@ Bool_t InitialStateHydjet::MultIni() {
   // calculation of  multiplicities of different particle species
   // according to the grand canonical approach
   GrandCanonical gc(15, fParams.fT, fParams.fMuB, fParams.fMuS, fParams.fMuI3);
-  GrandCanonical gc_ch(15, fParams.fT, fParams.fMuB, fParams.fMuS, fParams.fMuI3);
-  GrandCanonical gc_pi_th(15, fParams.fThFO, 0., 0., fParams.fMu_th_pip);
-  GrandCanonical gc_th_0(15, fParams.fThFO, 0., 0., 0.);
-  
-  // std::ofstream outMult("densities.txt");
-  // outMult<<"encoding    particle density      chemical potential "<<std::endl;
-  
+  GrandCanonical gcCh(15, fParams.fT, fParams.fMuB, fParams.fMuS, fParams.fMuI3);
+  GrandCanonical gcPiTh(15, fParams.fThFO, 0., 0., fParams.fMu_th_pip);
+  GrandCanonical gcTh0(15, fParams.fThFO, 0., 0., 0.);
   
   //effective volume for central     
   double dYl= 2 * fParams.fYlmax; //uniform distr. [-Ylmax; Ylmax]  
   if (fParams.fEtaType >0) dYl = TMath::Sqrt(2 * TMath::Pi()) * fParams.fYlmax ;  //Gaussian distr.                                                                            
   fVolEff = 2 * TMath::Pi() * fParams.fTau * dYl * (fParams.fR * fParams.fR)/TMath::Power((fParams.fUmax),2) * 
     ((fParams.fUmax)*TMath::SinH((fParams.fUmax))-TMath::CosH((fParams.fUmax))+ 1);
-  //  std::cout << "pion pointer = " << fDatabase->GetPDGParticle(211) << std::endl;
-  //  std::cout << "pion mass = " << (fDatabase->GetPDGParticle(211))->GetMass() << std::endl;
-  //  std::cout <<"central Effective volume = " << fVolEff << " [fm^3]" << std::endl;
   
-  Double_t particleDensity_pi_ch=0;
-  Double_t particleDensity_pi_th=0;
-  //  Double_t particleDensity_th_0=0;
+  Double_t particleDensityPiCh=0;
+  Double_t particleDensityPiTh=0;
   
   if(fParams.fThFO != fParams.fT && fParams.fThFO > 0){
-    //GrandCanonical gc_ch(15, fParams.fT, fParams.fMuB, fParams.fMuS, fParams.fMuI3);
-    //GrandCanonical gc_pi_th(15, fParams.fThFO, 0., 0., fParams.fMu_th_pip);
-    //GrandCanonical gc_th_0(15, fParams.fThFO, 0., 0., 0.);
-    particleDensity_pi_ch = gc_ch.ParticleNumberDensity(fDatabase->GetPDGParticle(211));
-    particleDensity_pi_th = gc_pi_th.ParticleNumberDensity(fDatabase->GetPDGParticle(211));
+    particleDensityPiCh = gcCh.ParticleNumberDensity(fDatabase->GetPDGParticle(211));
+    particleDensityPiTh = gcPiTh.ParticleNumberDensity(fDatabase->GetPDGParticle(211));
   }
 
   for(Int_t particleIndex = 0; particleIndex < fDatabase->GetNParticles(); particleIndex++) {
@@ -541,16 +491,12 @@ Bool_t InitialStateHydjet::MultIni() {
     
     //thermal f.o.
     if(fParams.fThFO != fParams.fT && fParams.fThFO > 0){
-      Double_t particleDensity_ch = gc_ch.ParticleNumberDensity(currParticle);
-      Double_t particleDensity_th_0 = gc_th_0.ParticleNumberDensity(currParticle);
-      Double_t numb_dens_bolt = particleDensity_pi_th*particleDensity_ch/particleDensity_pi_ch;               
-      //      Double_t coeff = ((currParticle->GetSpin() + 1.) * (currParticle->GetMass()) * 
-      //                       (currParticle->GetMass()) * fParams.fThFO / hbarc / hbarc / hbarc) /
-      //       (2. * TMath::Pi() * TMath::Pi()) * TMath::BesselK(2,(currParticle->GetMass())/fParams.fThFO);
-      //Double_t mumy = fParams.fThFO*TMath::Log(numb_dens_bolt/coeff);            
-      mu = fParams.fThFO*TMath::Log(numb_dens_bolt/particleDensity_th_0);
+      Double_t particleDensityCh = gcCh.ParticleNumberDensity(currParticle);
+      Double_t particleDensityTh0 = gcTh0.ParticleNumberDensity(currParticle);
+      Double_t numbDensBolt = particleDensityPiTh*particleDensityCh/particleDensityPiCh;               
+      mu = fParams.fThFO*TMath::Log(numbDensBolt/particleDensityTh0);
       if(abs(encoding)==211 || encoding==111)mu= fParams.fMu_th_pip; 
-      particleDensity = numb_dens_bolt;         
+      particleDensity = numbDensBolt;         
     }
     
     // set particle number density to zero for some species
@@ -580,14 +526,12 @@ Double_t InitialStateHydjet::SimpsonIntegrator2(Double_t a, Double_t b) {
   Int_t nsubIntervals=10000;
   Double_t h = (b - a)/nsubIntervals; //0-pi, phi
   Double_t s=0;
-  //  Double_t h2 = (fParams.fR)/nsubIntervals; //0-R maximal rB ?
-
   Double_t x = 0; //phi
   for(Int_t j = 1; j < nsubIntervals; j++) {
     x += h; // phi
     Double_t e = fParams.fEpsilon;
-    Double_t RsB = fParams.fR; //test: podstavit' *coefff_RB
-    Double_t rB = RsB *(TMath::Sqrt(1-e*e)/TMath::Sqrt(1+e*TMath::Cos(2*x))); //f-la7 rB    
+    Double_t rSB = fParams.fR; //test: podstavit' *coefff_RB
+    Double_t rB = rSB *(TMath::Sqrt(1-e*e)/TMath::Sqrt(1+e*TMath::Cos(2*x))); //f-la7 rB    
     Double_t sr = SimpsonIntegrator(0,rB,x);
     s += sr;
   }
@@ -617,8 +561,8 @@ Double_t InitialStateHydjet::SimpsonIntegrator(Double_t a, Double_t b, Double_t
 //f2=f(phi,r)
 Double_t InitialStateHydjet::f2(Double_t x, Double_t y) {
   // formula
-  Double_t RsB = fParams.fR; //test: podstavit' *coefff_RB
-  Double_t rhou =  fParams.fUmax * y / RsB;
+  Double_t rSB = fParams.fR; //test: podstavit' *coefff_RB
+  Double_t rhou =  fParams.fUmax * y / rSB;
   Double_t ff = y*TMath::CosH(rhou)*
     TMath::Sqrt(1+fParams.fDelta*TMath::Cos(2*x)*TMath::TanH(rhou)*TMath::TanH(rhou));
 //n_mu u^mu f-la 20
@@ -632,20 +576,14 @@ Double_t InitialStateHydjet::MidpointIntegrator2(Double_t a, Double_t b) {
   Int_t nsubIntervals2=1; 
   Double_t h = (b - a)/nsubIntervals; //0-pi , phi
   Double_t h2 = (fParams.fR)/nsubIntervals; //0-R maximal rB ?
-
   Double_t x = a + 0.5*h;
   Double_t y = 0;
-      
   Double_t t = f2(x,y);                    
   Double_t e = fParams.fEpsilon;
-
   for(Int_t j = 1; j < nsubIntervals; j++) {
     x += h; // integr  phi
-
-    Double_t RsB = fParams.fR; //test: podstavit' *coefff_RB
-    Double_t  rB = RsB *(TMath::Sqrt(1-e*e)/TMath::Sqrt(1+e*TMath::Cos(2*x))); //f-la7 rB
-
+    Double_t rSB = fParams.fR; //test: podstavit' *coefff_RB
+    Double_t  rB = rSB *(TMath::Sqrt(1-e*e)/TMath::Sqrt(1+e*TMath::Cos(2*x))); //f-la7 rB
     nsubIntervals2 = Int_t(rB / h2)+1;
     // integr R 
     y=0;
index 897ab7e..520e17a 100644 (file)
@@ -1,3 +1,10 @@
+// InitialStateHydjet is the class which controls the entire 
+// event simulation (input parameters, physics, output)
+// InitialStateHydjet inherits from class InitialState
+// MultIni() member function initializes PYQUEN and calculates the average soft multiplicities 
+// Initialize() member function calls the PYQUEN soubroutine for the hard part of the event
+// The resonance decay is performed by the function InitialState::Evolve()
+
 /*                                                                            
                                                                             
         Nikolai Amelin, Ludmila Malinina, Timur Pocheptsov (C) JINR/Dubna
 
 */
 
-#ifndef INITIALSTATEHYDJET_INCLUDED
-#define INITIALSTATEHYDJET_INCLUDED
+#ifndef INITIALSTATEHYDJET_H
+#define INITIALSTATEHYDJET_H
 
-#ifndef DATABASE_PDG
-#include "DatabasePDG.h"
-#endif
-#ifndef PARTICLE_INCLUDED
-#include "Particle.h"
-#endif
+//#ifndef DATABASE_PDG
+//#include "DatabasePDG.h"
+//#endif
+//#ifndef PARTICLE_INCLUDED
+//#include "Particle.h"
+//#endif
 #ifndef INITIAL_STATE
 #include "InitialState.h"
 #endif
 
+class ParticleAllocator;
+
 struct InitialParamsHydjet_t {
 
   Int_t fNevnt; //number of events
@@ -107,12 +116,18 @@ struct InitialParamsHydjet_t {
 
 class InitialStateHydjet : public InitialState {
  public:
+  InitialParamsHydjet_t fParams;             // the list of initial state parameters
+ private:
+  Double_t fVolEff;                           // the effective volume
+
+ public:
   InitialStateHydjet() : fParams(), fVolEff(0){};
   ~InitialStateHydjet() {};
   
   void SetVolEff(Double_t value) {fVolEff = value;}
-  Double_t GetVolEff() {return fVolEff;}
-  virtual Double_t GetTime() {return fParams.fDecay;}
+  Double_t GetVolEff() const {return fVolEff;}
+  //  virtual Double_t GetTime() {return fParams.fDecay;}
+  virtual Bool_t RunDecays() {return (fParams.fDecay>0 ? kTRUE : kFALSE);}
   virtual Int_t GetNev() {return fParams.fNevnt;}
   virtual Double_t GetWeakDecayLimit() {return fParams.fWeakDecay;}  
   
@@ -121,18 +136,12 @@ class InitialStateHydjet : public InitialState {
   virtual Bool_t MultIni();
   Bool_t IniOfThFreezeoutParameters();
  
-  Double_t f(Double_t);
-  Double_t f2(Double_t, Double_t);
+  Double_t f(Double_t x);
+  Double_t f2(Double_t x, Double_t y);
    
-  Double_t SimpsonIntegrator(Double_t, Double_t, Double_t);
-  Double_t SimpsonIntegrator2(Double_t, Double_t);
-  Double_t MidpointIntegrator2(Double_t, Double_t);
- public:
-  InitialParamsHydjet_t fParams;             // the list of initial state parameters  
-  
- private:
-  Double_t fVolEff;                           // the effective volume
-
+  Double_t SimpsonIntegrator(Double_t a, Double_t b, Double_t phi);
+  Double_t SimpsonIntegrator2(Double_t a, Double_t b);
+  Double_t MidpointIntegrator2(Double_t a, Double_t b);
 };
 
 #endif
index 6ced553..a0c4670 100644 (file)
@@ -22,13 +22,17 @@ Particle::Particle(ParticlePDG *prop):
   fParticleProperties(prop),
   fLastInteractionTime(0.),
   fInteractionNumber(0),
+  fPythiaStatusCode(-1),
   fLastMotherPdg(0),
   fType(0),
   fIndex(-1),
   fMotherIndex(-1),
-  fNDaughters(0)
+  fNDaughters(0),
+  fFirstDaughterIndex(-1),
+  fLastDaughterIndex(-1),
+  fDecayed(kFALSE)
 {
-  for(Int_t i=0; i<3; i++) fDaughterIndex[i] = -1;
+
 }
 
 Particle::Particle(ParticlePDG *prop, const TLorentzVector &pos, 
@@ -40,13 +44,16 @@ Particle::Particle(ParticlePDG *prop, const TLorentzVector &pos,
   fParticleProperties(prop),
   fLastInteractionTime(lit),
   fInteractionNumber(lin),
+  fPythiaStatusCode(-1),
   fLastMotherPdg(0),
   fType(type),
   fIndex(-1),
   fMotherIndex(-1),
-  fNDaughters(0)
+  fNDaughters(0),
+  fFirstDaughterIndex(-1),
+  fLastDaughterIndex(-1),
+  fDecayed(kFALSE)
 {
-  for(Int_t i=0; i<3; i++) fDaughterIndex[i] = -1;
 }
 
 Particle::Particle(ParticlePDG *prop, const TLorentzVector &pos, const TLorentzVector &mom,
@@ -59,36 +66,38 @@ Particle::Particle(ParticlePDG *prop, const TLorentzVector &pos, const TLorentzV
   fParticleProperties(prop),
   fLastInteractionTime(t),
   fInteractionNumber(n),
+  fPythiaStatusCode(-1),
   fLastMotherPdg(motherPdg),
   fType(ty),
   fIndex(-1),
   fMotherIndex(-1),
-  fNDaughters(0)
+  fNDaughters(0),
+  fFirstDaughterIndex(-1),
+  fLastDaughterIndex(-1),
+  fDecayed(kFALSE)
 {
-  for(Int_t i=0; i<3; i++) fDaughterIndex[i] = -1;
 }
 
 Particle::Particle(const Particle& copy) :
   fPosition(copy.Pos()),
-    fMomentum(copy.Mom()),
-    fLastMotherDecayCoor(copy.GetLastMotherDecayCoor()),
-    fLastMotherDecayMom(copy.GetLastMotherDecayMom()),
-    fParticleProperties(copy.Def()),
-    fLastInteractionTime(copy.GetLastInterTime()),
-    fInteractionNumber(copy.GetLastInterNumber()),
-    fLastMotherPdg(copy.GetLastMotherPdg()),
-    fType(copy.GetType()),
-    fIndex(copy.GetIndex()),
-    fMotherIndex(copy.GetMother()),
-    fNDaughters(copy.GetNDaughters())
+  fMomentum(copy.Mom()),
+  fLastMotherDecayCoor(copy.GetLastMotherDecayCoor()),
+  fLastMotherDecayMom(copy.GetLastMotherDecayMom()),
+  fParticleProperties(copy.Def()),
+  fLastInteractionTime(copy.GetLastInterTime()),
+  fInteractionNumber(copy.GetLastInterNumber()),
+  fPythiaStatusCode(copy.GetPythiaStatusCode()),
+  fLastMotherPdg(copy.GetLastMotherPdg()),
+  fType(copy.GetType()),
+  fIndex(copy.GetIndex()),
+  fMotherIndex(copy.GetMother()),
+  fNDaughters(copy.GetNDaughters()),
+  fFirstDaughterIndex(copy.GetFirstDaughterIndex()),
+  fLastDaughterIndex(copy.GetLastDaughterIndex()),
+  fDecayed(copy.GetDecayed())
 {
-      for(Int_t i=0; i<fNDaughters; i++) {
-        fDaughterIndex[i] = copy.GetDaughter(i);
-      }
 }
 
-
-
 Particle & Particle::operator=(const Particle& /*copy*/) {
   return *this;
 }
index aa1b4b8..7dc2cd0 100644 (file)
 #include "ParticlePDG.h"
 #endif
 #include <iostream>
-
+using namespace std;
 //class ParticlePDG;
 
 class Particle {
  public:
-
   Particle(const TLorentzVector &, const TLorentzVector &);
   Particle(const Particle& copy);
   Particle& operator=(const Particle&);
-
   virtual ~Particle() {};
   Particle(ParticlePDG *pdg = 0);
   Particle(ParticlePDG *pdg, const TLorentzVector &pos, const TLorentzVector &mom,
@@ -60,6 +57,9 @@ class Particle {
   TLorentzVector &Mom(){return fMomentum;}
   const TLorentzVector &Mom()const{return fMomentum;}
   TLorentzVector &Mom(const TLorentzVector &val){return fMomentum = val;}
+
+  void SetDecayed() {fDecayed = kTRUE;}
+  Bool_t GetDecayed() const {return fDecayed;}
                
   void Boost(const TVector3 &val){fMomentum.Boost(val);}
   void Boost(const TLorentzVector &val){fMomentum.Boost(val.BoostVector());}
@@ -85,33 +85,34 @@ class Particle {
 
   // aic(2008/08/08): functions added in order to enable tracking of mother/daughter particles by a unique index
   // The index coincides with the position of the particle in the secondaries list.
-  Int_t SetIndex() {fIndex = ++fLastIndex; return fIndex;}
+  Int_t SetIndex() {
+    fIndex = ++fLastIndex; 
+    return fIndex;
+  }
   Int_t GetIndex() const {return fIndex;}
   static Int_t GetLastIndex() {return fLastIndex;}
-  void InitIndexing() {fLastIndex = -1;}
+  static void InitIndexing() {
+    fLastIndex = -1;
+  }
   void SetMother(Int_t value) {fMotherIndex = value;}
   Int_t GetMother() const {return fMotherIndex;}
-  void SetDaughter(Int_t value) {
-    if(fNDaughters==3) {
-      std::cout << "Warning in Particle::SetDaughter() Already 3 daughters are set!! Check it out!!" << std::endl;
-      return;
-    }
-    fDaughterIndex[fNDaughters++] = value;
-  }
-  Int_t GetNDaughters() const {return fNDaughters;}
-  Int_t GetDaughter(Int_t value) const {
-    if(value<0 || value>fNDaughters-1) {
-      std::cout << "Warning in Particle::GetDaughter(Int_t) This particle has " << fNDaughters 
-               << " daughters. The argument must range from 0 to " << fNDaughters-1 << std::endl;
-      return -1;
-    }
-    return fDaughterIndex[value];
+  void SetFirstDaughterIndex(Int_t index) {fFirstDaughterIndex = index;}
+  void SetLastDaughterIndex(Int_t index) {fLastDaughterIndex = index;}
+  void SetPythiaStatusCode(Int_t code) {fPythiaStatusCode = code;}
+  Int_t GetPythiaStatusCode() const {return fPythiaStatusCode;}
+
+  //  Int_t GetNDaughters() const {return fNDaughters;}
+  Int_t GetNDaughters() const {
+    if(fFirstDaughterIndex==-1 || fLastDaughterIndex==-1) 
+      return 0;
+    else
+      return fLastDaughterIndex-fFirstDaughterIndex+1;
   }
+  Int_t GetFirstDaughterIndex() const {return fFirstDaughterIndex;}
+  Int_t GetLastDaughterIndex() const {return fLastDaughterIndex;}
 
-  // void SetLastMotherDecayCoor(TLorentzVector fLastMotherDecayCoor);
   TLorentzVector &SetLastMotherDecayCoor(const TLorentzVector &val){return fLastMotherDecayCoor = val;}
   const TLorentzVector &GetLastMotherDecayCoor()const{return fLastMotherDecayCoor;}
-  //  void SetLastMotherDecayMom(TLorentzVector fLastMotherDecayMom);
   TLorentzVector &SetLastMotherDecayMom(const TLorentzVector &val){return fLastMotherDecayMom = val;}
   const TLorentzVector &GetLastMotherDecayMom()const{return fLastMotherDecayMom;}
 
@@ -132,14 +133,16 @@ class Particle {
   ParticlePDG     *fParticleProperties;
   Double_t         fLastInteractionTime;
   Int_t            fInteractionNumber;
+  Int_t            fPythiaStatusCode;
   Int_t            fLastMotherPdg;
   Int_t            fType; //0-hydro, 1-jets
   Int_t            fIndex;                    // index (0 based) of particle in the final particle list which will contain both primaries and secondaries
   Int_t            fMotherIndex;              // index of the mother (-1 if its a primary particle)
   Int_t            fNDaughters;               // number of daughter particles (0 if the particle had not decayed)
-  Int_t            fDaughterIndex[3];         // array of indexes for the daughter particles (the indexes are -1 for non-existing daughters)
+  Int_t            fFirstDaughterIndex;       // index for the first daughter particle (-1 if non-existing)
+  Int_t            fLastDaughterIndex;        // index for the last daughter particle (-1 if non-existing)
+  Bool_t           fDecayed;                  // true if the decay procedure already applied
   static Int_t     fLastIndex;                // the last index assigned
-
 };
 
 Double_t S(const TLorentzVector &, const TLorentzVector &);
@@ -150,7 +153,7 @@ typedef std::list<Particle>::iterator LPIT_t;
 
 class ParticleAllocator {
  public:
-  ParticleAllocator() : fFreeNodes(0) {};
+  ParticleAllocator() {};
   void AddParticle(const Particle & particle, List_t & list);
   void FreeListNode(List_t & list, LPIT_t it);
   void FreeList(List_t & list);
index 39efae6..cf93ea9 100644 (file)
 #include "DecayChannel.h"
 #endif
 
-const Int_t kMaxDecayChannels = 20;
+const Int_t kMaxDecayChannels = 100;
 
 class ParticlePDG {
  public:
   ParticlePDG();
   ParticlePDG(Char_t* name, Int_t pdg, Double_t mass, Double_t width);
-
-  
   ~ParticlePDG();
   
   void AddChannel(DecayChannel* channel);
index 1661478..6b5bf68 100644 (file)
@@ -22,7 +22,7 @@
 #ifndef PARTICLE_INCLUDED
 #include "Particle.h"
 #endif
-#ifndef DATABASE_PDG
+#ifndef DATABASEPDG_H
 #include "DatabasePDG.h"
 #endif
 #ifndef PARTICLE_PDG
index 1e6c023..ad644e2 100644 (file)
@@ -18,7 +18,7 @@
 #ifndef NAEquationSolver_h
 #include "EquationSolver.h"
 #endif
-#ifndef DATABASE_PDG
+#ifndef DATABASEPDG_H
 #include "DatabasePDG.h"
 #endif