Skip to content
Snippets Groups Projects

Draft: Update Sexaquark anhillation simulation

Open Wouter Morren requested to merge wmorren/gauss-fork:jzhuo_sexaquark into jzhuo_sexaquark
2 files
+ 140
38
Compare changes
  • Side-by-side
  • Inline
Files
2
@@ -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);
Loading