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
+ 53
6
Compare changes
  • Side-by-side
  • Inline
Files
2
@@ -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
{
Loading