+Bool_t AliStack::IsStable(Int_t pdg) const
+{
+ //
+ // Decide whether particle (pdg) is stable
+ //
+
+
+ // All ions/nucleons are considered as stable
+ // Nuclear code is 10LZZZAAAI
+ if(pdg>1000000000)return kTRUE;
+
+ const Int_t kNstable = 18;
+ Int_t i;
+
+ Int_t pdgStable[kNstable] = {
+ kGamma, // Photon
+ kElectron, // Electron
+ kMuonPlus, // Muon
+ kPiPlus, // Pion
+ kKPlus, // Kaon
+ kK0Short, // K0s
+ kK0Long, // K0l
+ kProton, // Proton
+ kNeutron, // Neutron
+ kLambda0, // Lambda_0
+ kSigmaMinus, // Sigma Minus
+ kSigmaPlus, // Sigma Plus
+ 3312, // Xsi Minus
+ 3322, // Xsi
+ 3334, // Omega
+ kNuE, // Electron Neutrino
+ kNuMu, // Muon Neutrino
+ kNuTau // Tau Neutrino
+ };
+
+ 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 > 1) 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
+//
+
+ Int_t imo = p->GetFirstMother();
+ TParticle* pm = Particle(imo);
+ Int_t mpdg = TMath::Abs(pm->GetPdgCode());
+// Check for Sigma0
+ if ((mpdg == 3212) && (imo < GetNprimary())) return kTRUE;
+//
+// Check if it comes from a pi0 decay
+//
+// What about the pi0 Dalitz ??
+// if ((mpdg == kPi0) && (imo < GetNprimary())) return kTRUE;
+
+// Check if this is a heavy flavor decay product
+ 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 = pm->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 ?
+}