Geometry builder classes moved from base to sim.

[u/mrichter/AliRoot.git] / STEER / AliStack.cxx

diff --git a/STEER/AliStack.cxx b/STEER/AliStack.cxx
index 0295cba5b0071e59f030f4b3f3d6736356069f3f..8cea40392a617479c27f11944ea3161b64dc7157 100644 (file)
--- a/STEER/AliStack.cxx
+++ b/STEER/AliStack.cxx
@@ -29,6 +29,7 @@
 #include <TObjArray.h>
 #include <TParticle.h>
 #include <TParticlePDG.h>
+#include <TPDGCode.h>
 #include <TTree.h>
 
 #include "AliLog.h"
@@ -308,8 +309,7 @@ void AliStack::PurifyKine()
   if(fHgwmk+1 == fNtrack) return;
 
   // First pass, invalid Daughter information
-
-  for(i=0; i<fNtrack; i++) {
+   for(i=0; i<fNtrack; i++) {
     // Preset map, to be removed later
     if(i<=fHgwmk) map[i]=i ; 
     else {
@@ -368,6 +368,7 @@ void AliStack::PurifyKine()
       }
     }
   }
+
   
   // Now loop on all registered hit lists
   TList* hitLists = gAlice->GetMCApp()->GetHitLists();
@@ -377,6 +378,7 @@ void AliStack::PurifyKine()
   while((hitList = dynamic_cast<TCollection*>(next()))) {
     TIter nexthit(hitList);
     AliHit *hit;
+    
     while((hit = dynamic_cast<AliHit*>(nexthit()))) {
        hit->SetTrack(map[hit->GetTrack()]);
     }
@@ -431,7 +433,7 @@ void AliStack::ReorderKine()
   //
   // Howmany secondaries have been produced ?
   Int_t nNew = fNtrack - fHgwmk - 1;
-  
+    
   if (nNew > 0) {
       Int_t i, j;
       TObjArray &particles = *fParticleMap;
@@ -442,13 +444,7 @@ void AliStack::ReorderKine()
       
       for (i = 0; i < nNew; i++) {
          if (particles.At(fHgwmk + 1 + i)) {
-//           tmp[i] = (TParticle*) (particles.At(fHgwmk + 1 + i))->Clone();
              tmp[i] = (TParticle*) (particles.At(fHgwmk + 1 + i));
-
-//           if (((TParticle*) (particles.At(fHgwmk + 1 + i)))->TestBit(kKeepBit))
-//               tmp[i]->SetBit(kKeepBit);
-//           if (((TParticle*) (particles.At(fHgwmk + 1 + i)))->TestBit(kDoneBit))
-//               tmp[i]->SetBit(kDoneBit);
          } else {
              tmp[i] = 0x0;
          }
@@ -465,7 +461,7 @@ void AliStack::ReorderKine()
       // The outer loop is over parents, the inner over children.
       // -1 refers to the primary particle
       //
-      for (i = -1; i < nNew - 1; i++) {
+      for (i = -1; i < nNew-1; i++) {
          Int_t ipa;
          TParticle* parP;
          if (i == -1) {
@@ -480,9 +476,12 @@ void AliStack::ReorderKine()
              parP = tmp[i];
          }
           // Reset daughter information
+
+         Int_t idaumin = parP->GetFirstDaughter() - fHgwmk - 1;
+         Int_t idaumax = parP->GetLastDaughter()  - fHgwmk - 1;
          parP->SetFirstDaughter(-1);
          parP->SetLastDaughter(-1);
-         for (j = i + 1; j < nNew; j++) {
+         for (j = idaumin; j <= idaumax; j++) {
               // Skip deleted particles
              if (!tmp[j])        continue;
               // Skip particles already handled
@@ -1045,3 +1044,100 @@ void AliStack::SetEventFolderName(const char* foldname)
  //Sets event folder name
  fEventFolderName = foldname;
 }
+
+Bool_t AliStack::IsStable(Int_t pdg) const
+{
+//
+// Decide whether particle (pdg) is stable
+//
+
+    const Int_t kNstable = 14;
+    Int_t i;
+
+    Int_t pdgStable[kNstable] = {
+       kGamma,             // Photon
+       kElectron,          // Electron
+       kMuonPlus,          // Muon 
+       kPiPlus,            // Pion
+       kKPlus,             // Kaon
+       kProton,            // Proton 
+       kNeutron,           // Neutron
+       kLambda0,           // Lambda_0
+       kSigmaMinus,        // Sigma Minus
+       kSigma0,            // Sigma_0
+       kSigmaPlus,         // Sigma Plus
+       3312,               // Xsi Minus 
+       3322,               // Xsi 
+       3334,               // Omega
+    };
+    
+    Bool_t isStable = kFALSE;
+    for (i = 0; i < kNstable; i++) {
+       if (pdg == TMath::Abs(pdgStable[i])) {
+           isStable = kTRUE;
+           break;
+       }
+    }
+
+    return isStable;
+}
+
+Bool_t AliStack::IsPhysicalPrimary(Int_t index)
+{
+    //
+    // Test if a particle is a physical primary according to the following definition:
+    // Particles produced in the collision including products of strong and
+    // electromagnetic decay and excluding feed-down from weak decays of strange
+    // particles.
+    //
+    TParticle* p = Particle(index);
+    Int_t ist = p->GetStatusCode();
+    
+    //
+    // Initial state particle
+    if (ist > 20) return kFALSE;
+    
+    Int_t pdg = TMath::Abs(p->GetPdgCode());
+    
+    if (!IsStable(pdg)) return kFALSE;
+    
+    if (index < GetNprimary()) {
+//
+// Particle produced by generator
+       return kTRUE;
+    } else {
+//
+// Particle produced during transport
+//
+// Check if this is a heavy flavor decay product
+       Int_t imo =  p->GetFirstMother();
+       TParticle* pm  = Particle(imo);
+       Int_t mpdg = TMath::Abs(pm->GetPdgCode());
+       Int_t mfl  = Int_t (mpdg / TMath::Power(10, Int_t(TMath::Log10(mpdg))));
+       //
+       // Light hadron
+       if (mfl < 4) return kFALSE;
+       
+       //
+       // Heavy flavor hadron produced by generator
+       if (imo <  GetNprimary()) {
+           return kTRUE;
+       }
+       
+       // To be sure that heavy flavor has not been produced in a secondary interaction
+       // Loop back to the generated mother
+       while (imo >=  GetNprimary()) {
+           imo = p->GetFirstMother();
+           pm  =  Particle(imo);
+       }
+       mpdg = TMath::Abs(pm->GetPdgCode());
+       mfl  = Int_t (mpdg / TMath::Power(10, Int_t(TMath::Log10(mpdg))));
+
+       if (mfl < 4) {
+           return kFALSE;
+       } else {
+           return kTRUE;
+       } 
+    } // produced by generator ?
+} 
+