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