From 91fe033caac38390d51b0bcc27138dcfbdf24a6f Mon Sep 17 00:00:00 2001
From: Wouter Max Morren <wmorren@cern.ch>
Date: Wed, 28 Feb 2024 14:41:28 +0100
Subject: [PATCH 1/3] Initial commit

---
 .../src/G4AntiSexaquarkAnnihilation.hpp       | 57 +++++++++++++++++--
 Sim/GaussPhysics/src/G4Sexaquark.cpp          |  2 +-
 2 files changed, 53 insertions(+), 6 deletions(-)

diff --git a/Sim/GaussPhysics/src/G4AntiSexaquarkAnnihilation.hpp b/Sim/GaussPhysics/src/G4AntiSexaquarkAnnihilation.hpp
index 013d53e32..4a023e11a 100644
--- a/Sim/GaussPhysics/src/G4AntiSexaquarkAnnihilation.hpp
+++ b/Sim/GaussPhysics/src/G4AntiSexaquarkAnnihilation.hpp
@@ -108,7 +108,7 @@ namespace
   template<typename F>
   std::optional<const G4Element*> GetInteractionElement(const G4Step& aStep, F const& f)
   {
-    const auto material = aStep.GetPreStepPoint()->GetMaterial();
+        const auto material = aStep.GetPreStepPoint()->GetMaterial();
     const auto elements = material->GetElementVector();
 
     for (unsigned i=0; i<elements->size(); i++)
@@ -139,14 +139,61 @@ G4bool G4AntiSexaquarkAnnihilation_t<InteractionType, InteractionConfig>::IsAppl
 template<G4AntiSexaquarkAnnihilationType InteractionType, G4AntiSexaquarkAnnihilationConfig InteractionConfig>
 G4double G4AntiSexaquarkAnnihilation_t<InteractionType, InteractionConfig>::GetMeanFreePath( const G4Track& aTrack, G4double, G4ForceCondition*)
 {
+  //When I decrease the cross section from 1e-16 to 1e-18 the error occured, meaning that suddenly S is able to travel further, 
+  //where it encounters a material for which GetA is not working.
+
+  //Physical volume : /dd/Geometry/BeforeMagnetRegion/Rich1/lvRich1SubMaster#pvUX851InRich1SubMaster#pvUX851Cyl04 (/dd/Materials/Pipe/PipeBeTV56)
+  //Physical volume : /dd/Geometry/BeforeMagnetRegion/VP/lvVP#pvVPRight#pvModule22WithSupport#pvModule22#pvLadder22_2#pvChips22_2#pvGlue0 (/dd/Materials/VP/Stycast)
+
+  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 fNumberOfElements = 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 (fNumberOfElements > 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 (std::count(theAtomsVector, theAtomsVector + fNumberOfElements, 0) > 0 ){ // composition by fractional mass
+        const auto theFractionVector = material->GetFractionVector();
+
+        for (G4int i=0; i<fNumberOfElements; ++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<fNumberOfElements; ++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
   {
diff --git a/Sim/GaussPhysics/src/G4Sexaquark.cpp b/Sim/GaussPhysics/src/G4Sexaquark.cpp
index 9895e2ebc..c77344131 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") ; }
-- 
GitLab


From 4ccd941c9e20a8c29ff74d286a29ea4795937e46 Mon Sep 17 00:00:00 2001
From: Wouter Max Morren <wmorren@cern.ch>
Date: Tue, 12 Mar 2024 15:44:57 +0100
Subject: [PATCH 2/3] Randomly select daughter particle & disallow 3 body
 anhillation

---
 .../src/G4AntiSexaquarkAnnihilation.hpp       | 151 ++++++++++++------
 1 file changed, 103 insertions(+), 48 deletions(-)

diff --git a/Sim/GaussPhysics/src/G4AntiSexaquarkAnnihilation.hpp b/Sim/GaussPhysics/src/G4AntiSexaquarkAnnihilation.hpp
index 4a023e11a..1e751ac13 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"
@@ -106,17 +109,76 @@ public:
 namespace
 {
   template<typename F>
-  std::optional<const G4Element*> GetInteractionElement(const G4Step& aStep, F const& f)
+  std::optional<const G4Isotope*> GetInteractionElement(const G4Step& aStep, F const& f, const bool& Proton)
   {
-        const auto material = aStep.GetPreStepPoint()->GetMaterial();
+    /*
+    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) { pAtom[0]=1; } // No composition
+    else {
+      double totAtom = 0; // Normalisation factor
+      for (unsigned i=0; i<nElements; i++)
+      {
+        double Z = (*elements)[i]->GetZ();
+        double N = (*elements)[i]->GetN();
+        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 (f(*element) && 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 Z = (*theIsotopeVector)[i]->GetZ();
+            double N = (*theIsotopeVector)[i]->GetN();
+            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,47 +201,39 @@ G4bool G4AntiSexaquarkAnnihilation_t<InteractionType, InteractionConfig>::IsAppl
 template<G4AntiSexaquarkAnnihilationType InteractionType, G4AntiSexaquarkAnnihilationConfig InteractionConfig>
 G4double G4AntiSexaquarkAnnihilation_t<InteractionType, InteractionConfig>::GetMeanFreePath( const G4Track& aTrack, G4double, G4ForceCondition*)
 {
-  //When I decrease the cross section from 1e-16 to 1e-18 the error occured, meaning that suddenly S is able to travel further, 
-  //where it encounters a material for which GetA is not working.
-
-  //Physical volume : /dd/Geometry/BeforeMagnetRegion/Rich1/lvRich1SubMaster#pvUX851InRich1SubMaster#pvUX851Cyl04 (/dd/Materials/Pipe/PipeBeTV56)
-  //Physical volume : /dd/Geometry/BeforeMagnetRegion/VP/lvVP#pvVPRight#pvModule22WithSupport#pvModule22#pvLadder22_2#pvChips22_2#pvGlue0 (/dd/Materials/VP/Stycast)
-
   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 (!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 theElementVector = material->GetElementVector(); //vector of pointers to elements constituing this material
-    const G4int fNumberOfElements = material->GetNumberOfElements();
+    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 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
 
-    if (fNumberOfElements > 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 (std::count(theAtomsVector, theAtomsVector + fNumberOfElements, 0) > 0 ){ // composition by fractional mass
+      // 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<fNumberOfElements; ++i) {
+        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<fNumberOfElements; ++i) {
+      else { // Composition by atom count
+        for (G4int i=0; i<nElements; ++i) {
           N += (*theElementVector)[i]->GetN() * theAtomsVector[i];
           Z += (*theElementVector)[i]->GetZ() * theAtomsVector[i];
         }
@@ -187,7 +241,7 @@ G4double G4AntiSexaquarkAnnihilation_t<InteractionType, InteractionConfig>::GetM
         else if ( Neutron ) { f = (N-Z)/N; }
       }
     }
-    else {// no composition
+    else { // No composition
       N = (*theElementVector)[0]->GetN(); 
       Z = material->GetZ();
       if (Proton)         { f = Z/N; }
@@ -224,8 +278,8 @@ 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);
+    const auto HasProton = [](const G4Element & element){ return (element.GetZ() >= 1.); }; // Isn't 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();
@@ -233,8 +287,8 @@ G4VParticleChange* G4AntiSexaquarkAnnihilation_t<InteractionType, InteractionCon
   }
   else if constexpr ( InteractionType == G4AntiSexaquarkAnnihilationType::NeutronInteraction )
   {
-    const auto HasNeutron = [](const G4Element & element){ return (element.GetN() >= 1.); };
-    const auto InteractionElement = GetInteractionElement(aStep, HasNeutron);
+    const auto HasNeutron = [](const G4Element & element){ return (element.GetN() - element.GetZ() >= 1.); }; // Here is isn't always the case indeed
+    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();
@@ -242,27 +296,28 @@ G4VParticleChange* G4AntiSexaquarkAnnihilation_t<InteractionType, InteractionCon
   }
   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() > 1.) && (element.GetN() - element.GetZ() > 0.); };
+    const auto InteractionElement = GetInteractionElement(aStep, ValidAtom, true); //The first atom of the material is now always chosen.
     if ( !InteractionElement.has_value() ) { return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep ); }
-    const auto TargetA = InteractionElement.value()->GetA();
-    const auto TargetZ = InteractionElement.value()->GetZ();
+    const auto Z = InteractionElement.value()->GetZ();
+    const auto N = InteractionElement.value()->GetN(); 
+    fmt::print("Z:{}, N:{}\n", Z, N);
     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);
+    const auto ValidAtom = [](const G4Element & element){ return (element.GetZ() > 0.) && (element.GetN() - element.GetZ() > 1.); };
+    const auto InteractionElement = GetInteractionElement(aStep, ValidAtom, false);
     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] = 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
   {
@@ -275,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);
 
-- 
GitLab


From 021efb3b4c3fe3bc713a66bdf3d03e319d969237 Mon Sep 17 00:00:00 2001
From: Wouter Max Morren <wmorren@cern.ch>
Date: Wed, 13 Mar 2024 15:10:52 +0100
Subject: [PATCH 3/3] Some simplification, but left some changes as comments
 for now

---
 .../src/G4AntiSexaquarkAnnihilation.hpp       | 66 +++++++++----------
 1 file changed, 33 insertions(+), 33 deletions(-)

diff --git a/Sim/GaussPhysics/src/G4AntiSexaquarkAnnihilation.hpp b/Sim/GaussPhysics/src/G4AntiSexaquarkAnnihilation.hpp
index 1e751ac13..8154fb77b 100644
--- a/Sim/GaussPhysics/src/G4AntiSexaquarkAnnihilation.hpp
+++ b/Sim/GaussPhysics/src/G4AntiSexaquarkAnnihilation.hpp
@@ -104,12 +104,12 @@ public:
 };
 
 ///
-/// Implementation
+/// Implementation for charged and neutral interaction only
 ///
-namespace
-{
-  template<typename F>
-  std::optional<const G4Isotope*> GetInteractionElement(const G4Step& aStep, F const& f, const bool& Proton)
+//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.
@@ -118,29 +118,30 @@ namespace
     */
     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) { pAtom[0]=1; } // No composition
-    else {
+    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 Z = (*elements)[i]->GetZ();
         double N = (*elements)[i]->GetN();
-        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);}
+        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];
         }
-        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
@@ -151,7 +152,7 @@ namespace
     for (unsigned i=0; i<nElements; i++)
     {
       const auto element = elements->at(i);
-      if (f(*element) && r1 < pAtom[i]){
+      if (r1 < pAtom[i]){
         const auto theIsotopeVector = element->GetIsotopeVector();
         const auto theRelativeAbundanceVector = element->GetRelativeAbundanceVector();
         const size_t nIsotopes = element->GetNumberOfIsotopes();
@@ -162,8 +163,8 @@ namespace
         else {
           double totIsotope = 0; // Normalisation factor
           for (unsigned i=0; i<nIsotopes; i++) {
-            double Z = (*theIsotopeVector)[i]->GetZ();
             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];
@@ -183,7 +184,7 @@ namespace
     }
     return std::nullopt;
   }
-}
+//}
 
 
 //
@@ -278,30 +279,29 @@ G4VParticleChange* G4AntiSexaquarkAnnihilation_t<InteractionType, InteractionCon
 
   if constexpr ( InteractionType == G4AntiSexaquarkAnnihilationType::ProtonInteraction )
   {
-    const auto HasProton = [](const G4Element & element){ return (element.GetZ() >= 1.); }; // Isn't this always the case?
-    const auto InteractionElement = GetInteractionElement(aStep, HasProton, true);
-    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() - element.GetZ() >= 1.); }; // Here is isn't always the case indeed
-    const auto InteractionElement = GetInteractionElement(aStep, HasNeutron, false);
-    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() > 1.) && (element.GetN() - element.GetZ() > 0.); };
-    const auto InteractionElement = GetInteractionElement(aStep, ValidAtom, true); //The first atom of the material is now always chosen.
+    //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 Z = InteractionElement.value()->GetZ();
     const auto N = InteractionElement.value()->GetN(); 
-    fmt::print("Z:{}, N:{}\n", Z, N);
+    const auto Z = InteractionElement.value()->GetZ();
     const auto ionTable = G4IonTable::GetIonTable();
     daughters_definitions[0] = G4AntiXiMinus::Definition();
     daughters_definitions[1] = ionTable->GetIon(Z-1, N);
@@ -309,9 +309,9 @@ G4VParticleChange* G4AntiSexaquarkAnnihilation_t<InteractionType, InteractionCon
   }
   else if constexpr ( InteractionType == G4AntiSexaquarkAnnihilationType::NeutralAtomicInteraction )
   {
-    const auto ValidAtom = [](const G4Element & element){ return (element.GetZ() > 0.) && (element.GetN() - element.GetZ() > 1.); };
-    const auto InteractionElement = GetInteractionElement(aStep, ValidAtom, false);
-    if ( !InteractionElement.has_value() ) { return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep ); }
+    //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();
-- 
GitLab