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") ; }