diff --git a/Sim/GaussPhysics/src/G4AntiSexaquarkAnnihilation.hpp b/Sim/GaussPhysics/src/G4AntiSexaquarkAnnihilation.hpp
index 013d53e32c2c69c227d01d62f91eb6704000180c..8154fb77b14d4df78d2e77aa64e054e074a1c2da 100644
--- a/Sim/GaussPhysics/src/G4AntiSexaquarkAnnihilation.hpp
+++ b/Sim/GaussPhysics/src/G4AntiSexaquarkAnnihilation.hpp
@@ -10,6 +10,9 @@
 
 // Geant4
 #include "G4IonTable.hh"
+#include "G4VIsotopeTable.hh"
+#include "G4Isotope.hh"
+#include "G4IsotopeVector.hh"
 #include "G4ForceCondition.hh"
 #include "Randomize.hh"
 #include "G4Step.hh"
@@ -101,27 +104,87 @@ public:
 };
 
 ///
-/// Implementation
+/// Implementation for charged and neutral interaction only
 ///
-namespace
-{
-  template<typename F>
-  std::optional<const G4Element*> GetInteractionElement(const G4Step& aStep, F const& f)
+//namespace
+//{
+  //template<typename F>
+  std::optional<const G4Isotope*> GetInteractionElement(const G4Step& aStep, const bool& Proton) //F const& f,
   {
+    /*
+    Need to choose an atom and isotope as daughter particle for charged and neutral atomic interaction.
+    The probability of hitting a specific atom can be calculated using the (mean) number of protons or neutrons and the atom count or fraction. 
+    Then the probability of hitting a specific isotope is dependent on its abundance and number of protons or neutrons inside.
+    */
     const auto material = aStep.GetPreStepPoint()->GetMaterial();
     const auto elements = material->GetElementVector();
+    const size_t nElements = material->GetNumberOfElements();
+    const auto theFractionVector = material->GetFractionVector();
+    const auto theAtomsVector = material->GetAtomsVector(); // Vector of atom count of each element
+
+    // Calculate probablity to interact with specific atom
+    std::vector<double> pAtom(nElements); // Probability vector
+    if (nElements==1 && ((*elements)[0]->GetZ() >= 2. || !Proton)) { pAtom[0]=1; } // No composition & no Hydrogen charged interation
+    else if (nElements >1) {
+      double totAtom = 0; // Normalisation factor
+      for (unsigned i=0; i<nElements; i++)
+      {
+        double N = (*elements)[i]->GetN();
+        double Z = (*elements)[i]->GetZ();
+        if (Z >= 2. || !Proton) { // Disallow 3 body charged anhillation with Hydrogen 
+          if (std::count(theAtomsVector, theAtomsVector + nElements, 0) > 0 ){ // Composition by fractional mass
+            if (Proton) {pAtom[i] = theFractionVector[i] * Z;}
+            else {pAtom[i] = theFractionVector[i] * (N - Z);}
+          }
+          else{ // Composition by atom count or single atom
+            if (Proton) {pAtom[i] = theAtomsVector[i] * Z;}
+            else {pAtom[i] = theAtomsVector[i] * (N - Z);}
+          }
+          totAtom+=pAtom[i];
+        }
+      }
+      for (unsigned i=0; i<nElements; i++) { pAtom[i] /= totAtom; } // Normalise
+      for (unsigned i=1; i<nElements; i++) { pAtom[i] += pAtom[i-1]; } // Cumulate
+    }
 
-    for (unsigned i=0; i<elements->size(); i++)
+    // Choose atom
+    double r1 = G4UniformRand();
+    for (unsigned i=0; i<nElements; i++)
     {
       const auto element = elements->at(i);
-      if (f(*element))
-      {
-        return element;
+      if (r1 < pAtom[i]){
+        const auto theIsotopeVector = element->GetIsotopeVector();
+        const auto theRelativeAbundanceVector = element->GetRelativeAbundanceVector();
+        const size_t nIsotopes = element->GetNumberOfIsotopes();
+        
+        // Calculate probability to interact with specific isotope
+        std::vector<double> pIsotope(nIsotopes); // Probability vector
+        if (nIsotopes==1) { pIsotope[0]=1; }
+        else {
+          double totIsotope = 0; // Normalisation factor
+          for (unsigned i=0; i<nIsotopes; i++) {
+            double N = (*theIsotopeVector)[i]->GetN();
+            double Z = (*theIsotopeVector)[i]->GetZ();
+            if (Proton){pIsotope[i] = theRelativeAbundanceVector[i] * Z; }
+            else {pIsotope[i] = theRelativeAbundanceVector[i] * (N - Z); }
+            totIsotope += pIsotope[i];
+          }
+          for (unsigned i=0; i<nIsotopes; i++) { pIsotope[i] /= totIsotope; } // Normalise
+          for (unsigned i=1; i<nIsotopes; i++) { pIsotope[i] += pIsotope[i-1]; } // Cumulate
+        }
+
+        // Choose isotope
+        double r2 = G4UniformRand();
+        for (unsigned i=0; i<nIsotopes; i++){
+          if (r2 < pIsotope[i]){
+            return (*theIsotopeVector)[i];
+          }
+        }
       }
     }
     return std::nullopt;
   }
-}
+//}
 
 
 //
@@ -139,14 +202,53 @@ G4bool G4AntiSexaquarkAnnihilation_t<InteractionType, InteractionConfig>::IsAppl
 template<G4AntiSexaquarkAnnihilationType InteractionType, G4AntiSexaquarkAnnihilationConfig InteractionConfig>
 G4double G4AntiSexaquarkAnnihilation_t<InteractionType, InteractionConfig>::GetMeanFreePath( const G4Track& aTrack, G4double, G4ForceCondition*)
 {
+  const G4bool Proton = InteractionType == G4AntiSexaquarkAnnihilationType::ProtonInteraction || InteractionType == G4AntiSexaquarkAnnihilationType::ChargeAtomicInteraction;
+  const G4bool Neutron = InteractionType == G4AntiSexaquarkAnnihilationType::NeutronInteraction || InteractionType == G4AntiSexaquarkAnnihilationType::NeutralAtomicInteraction;
+  if (!Proton && !Neutron) { throw GiGaException ("G4AntiSexaquarkAnnihilation::PostStepDoIt: Invalid interaction type"); }
+
   if constexpr (InteractionConfig == G4AntiSexaquarkAnnihilationConfig::FromCrossSection)
   {
     const auto material = aTrack.GetMaterial();
     const auto density = material->GetDensity();
-    const auto A = material->GetA();
-    const auto n = density / A * Gaudi::Units::Avogadro;
-    const auto sigma = m_configuration_value;
-    return 1./(n*sigma);
+    const auto theElementVector = material->GetElementVector(); // Vector of pointers to elements constituing this material
+    const G4int nElements = material->GetNumberOfElements();
+    G4double N = 0.; // Mass number
+    G4double Z = 0.; // Proton number
+    G4double f = 0.; // Fraction of protons or neutrons in material
+    G4double n = density * Gaudi::Units::Avogadro / 1e24; // Number of nucleons per cm3 
+    const auto sigma = m_configuration_value * 1e24; // Cross-section
+    //1e24 to make sure n and sigma are somewhat same order of magnitude to reduce numerical errors during multiplication
+    
+    if (nElements > 1) { // Composite material
+      const auto theAtomsVector = material->GetAtomsVector(); // Vector of atom count of each element
+
+      // Sums all zero atom counts between first pointer and last pointer, if 0 in list use fraction vector
+      if (std::count(theAtomsVector, theAtomsVector + nElements, 0) > 0 ){ // Composition by fractional mass
+        const auto theFractionVector = material->GetFractionVector();
+
+        for (G4int i=0; i<nElements; ++i) {
+          N = (*theElementVector)[i]->GetN();
+          Z = (*theElementVector)[i]->GetZ();
+          if (Proton)       { f += theFractionVector[i] * Z / N; }
+          else if (Neutron) { f += theFractionVector[i] * (N - Z) / N; }
+        }
+      }
+      else { // Composition by atom count
+        for (G4int i=0; i<nElements; ++i) {
+          N += (*theElementVector)[i]->GetN() * theAtomsVector[i];
+          Z += (*theElementVector)[i]->GetZ() * theAtomsVector[i];
+        }
+        if (Proton)         { f = Z/N; }
+        else if ( Neutron ) { f = (N-Z)/N; }
+      }
+    }
+    else { // No composition
+      N = (*theElementVector)[0]->GetN(); 
+      Z = material->GetZ();
+      if (Proton)         { f = Z/N; }
+      else if ( Neutron ) { f = (N-Z)/N; }
+    } 
+    return 1./(n*sigma*f);
   }
   else
   {
@@ -177,45 +279,45 @@ G4VParticleChange* G4AntiSexaquarkAnnihilation_t<InteractionType, InteractionCon
 
   if constexpr ( InteractionType == G4AntiSexaquarkAnnihilationType::ProtonInteraction )
   {
-    const auto HasProton = [](const G4Element & element){ return (element.GetZ() >= 1.); };
-    const auto InteractionElement = GetInteractionElement(aStep, HasProton);
-    if ( !InteractionElement.has_value() ) { return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep ); }
+    //const auto HasProton = [](const G4Element & element){ return (element.GetZ() >= 1.); }; // This always the case
+    //const auto InteractionElement = GetInteractionElement(aStep, HasProton, true);
+    //if ( !InteractionElement.has_value() ) { return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep ); }
     daughters_definitions[0] = G4AntiLambda::Definition();
     daughters_definitions[1] = G4KaonPlus  ::Definition();
     target_mass = G4Proton::Definition()->GetPDGMass();
   }
   else if constexpr ( InteractionType == G4AntiSexaquarkAnnihilationType::NeutronInteraction )
   {
-    const auto HasNeutron = [](const G4Element & element){ return (element.GetN() >= 1.); };
-    const auto InteractionElement = GetInteractionElement(aStep, HasNeutron);
-    if ( !InteractionElement.has_value() ) { return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep ); }
+    //const auto HasNeutron = [](const G4Element & element){ return (element.GetN() - element.GetZ() >= 1.); }; 
+    //const auto InteractionElement = GetInteractionElement(aStep, HasNeutron, false);
+    //if ( !InteractionElement.has_value() ) { return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep ); }
     daughters_definitions[0] = G4AntiLambda   ::Definition();
     daughters_definitions[1] = G4KaonZeroShort::Definition();
     target_mass = G4Neutron::Definition()->GetPDGMass();
   }
   else if constexpr ( InteractionType == G4AntiSexaquarkAnnihilationType::ChargeAtomicInteraction )
   {
-    const auto ValidAtom = [](const G4Element & element){ return (element.GetZ() > 0.) && (element.GetN() > 0.); };
-    const auto InteractionElement = GetInteractionElement(aStep, ValidAtom);
+    //const auto ValidAtom = [](const G4Element & element){ return (element.GetZ() >= 2.) && (element.GetN() - element.GetZ() > 0.); }; // 3 body anhillation with H kinemtically not allowed (ignoring D)
+    const auto InteractionElement = GetInteractionElement(aStep, true); 
     if ( !InteractionElement.has_value() ) { return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep ); }
-    const auto TargetA = InteractionElement.value()->GetA();
-    const auto TargetZ = InteractionElement.value()->GetZ();
+    const auto N = InteractionElement.value()->GetN(); 
+    const auto Z = InteractionElement.value()->GetZ();
     const auto ionTable = G4IonTable::GetIonTable();
     daughters_definitions[0] = G4AntiXiMinus::Definition();
-    daughters_definitions[1] = ionTable->GetIon(TargetZ-1, TargetA);
-    target_mass = ionTable->GetIonMass(TargetZ, TargetA);
+    daughters_definitions[1] = ionTable->GetIon(Z-1, N);
+    target_mass = ionTable->GetIonMass(Z, N);
   }
   else if constexpr ( InteractionType == G4AntiSexaquarkAnnihilationType::NeutralAtomicInteraction )
   {
-    const auto ValidAtom = [](const G4Element & element){ return (element.GetZ() > 0.) && (element.GetN() > 0.); };
-    const auto InteractionElement = GetInteractionElement(aStep, ValidAtom);
-    if ( !InteractionElement.has_value() ) { return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep ); }
-    const auto TargetA = InteractionElement.value()->GetA();
-    const auto TargetZ = InteractionElement.value()->GetZ();
+    //const auto ValidAtom = [](const G4Element & element){ return (element.GetZ() >= 1.) && (element.GetN() - element.GetZ() >= 1.); }; 
+    const auto InteractionElement = GetInteractionElement(aStep, false);
+    //if ( !InteractionElement.has_value() ) { return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep ); }
+    const auto N = InteractionElement.value()->GetN();
+    const auto Z = InteractionElement.value()->GetZ();
     const auto ionTable = G4IonTable::GetIonTable();
     daughters_definitions[0] = G4AntiXiZero::Definition();
-    daughters_definitions[1] = ionTable->GetIon(TargetZ, TargetA-1);
-    target_mass = ionTable->GetIonMass(TargetZ, TargetA);
+    daughters_definitions[1] = ionTable->GetIon(Z-1, N);
+    target_mass = ionTable->GetIonMass(Z, N);
   }
   else
   {
@@ -228,14 +330,14 @@ G4VParticleChange* G4AntiSexaquarkAnnihilation_t<InteractionType, InteractionCon
 
   // Get daugther masses
   std::vector<G4double> daughter_masses =
-  {
-    daughters_definitions[0]->GetPDGMass(), daughters_definitions[1]->GetPDGMass()
-  };
+    {
+      daughters_definitions[0]->GetPDGMass(), daughters_definitions[1]->GetPDGMass()
+    };
 
   // Generate PHSP end-states
   std::vector<G4LorentzVector> end_states;
   m_PHSP_generator->Generate( init_state, daughter_masses, end_states );
-
+  
   // Stop and kill input track
   aParticleChange.ProposeTrackStatus(fStopAndKill);
 
diff --git a/Sim/GaussPhysics/src/G4Sexaquark.cpp b/Sim/GaussPhysics/src/G4Sexaquark.cpp
index 9895e2ebc843c583634ea81900ef4a38eb21e570..c77344131a8fce7c4dcc9317133ae16a22826f54 100644
--- a/Sim/GaussPhysics/src/G4Sexaquark.cpp
+++ b/Sim/GaussPhysics/src/G4Sexaquark.cpp
@@ -59,7 +59,7 @@ G4Sexaquark* G4Sexaquark::Definition(){ return s_instance; }
 G4Sexaquark* G4Sexaquark::Create(const double mass)
 {
   if ( 0 != s_instance ) { return s_instance ; }
-  // check the resence of the paricle in particle table:
+  // check the presence of the paricle in particle table:
   G4ParticleTable* table = G4ParticleTable::GetParticleTable();
   if ( 0 == table )
   { throw GiGaException ("G4Sexaquark::Create: Invalid G4ParticleTable") ; }