From 7ca2a7cd9746f01f6f7bbcfde3dbc4f257f124bc Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Sat, 22 Jun 2019 16:13:35 +0100
Subject: [PATCH] added medium composition to sibyll interaction constructor
---
Documentation/Examples/boundary_example.cc | 15 ++++++++++-
Documentation/Examples/cascade_example.cc | 16 +++++++++++-
Documentation/Examples/vertical_EAS.cc | 15 ++++++++++-
Processes/Sibyll/Interaction.cc | 30 ++++++----------------
Processes/Sibyll/Interaction.h | 11 ++++++--
Processes/Sibyll/testSibyll.cc | 13 ++++++----
6 files changed, 68 insertions(+), 32 deletions(-)
diff --git a/Documentation/Examples/boundary_example.cc b/Documentation/Examples/boundary_example.cc
index 7b990016..52ba2f38 100644
--- a/Documentation/Examples/boundary_example.cc
+++ b/Documentation/Examples/boundary_example.cc
@@ -115,11 +115,24 @@ int main() {
universe.AddChild(std::move(outerMedium));
+ auto const allElementsInUniverse = std::invoke([&]() {
+ std::set<particles::Code> allElementsInUniverse;
+ auto collectElements = [&](auto& vtn) {
+ if (vtn.HasModelProperties()) {
+ auto const& comp =
+ vtn.GetModelProperties().GetNuclearComposition().GetComponents();
+ for (auto const c : comp) allElementsInUniverse.insert(c);
+ }
+ };
+ universe.walk(collectElements);
+ return allElementsInUniverse;
+ });
+
// setup processes, decays and interactions
tracking_line::TrackingLine tracking;
random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
- process::sibyll::Interaction sibyll;
+ process::sibyll::Interaction sibyll(allElementsInUniverse);
process::sibyll::Decay decay;
process::particle_cut::ParticleCut cut(20_GeV);
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index c2cc3c4f..f07089a1 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -92,6 +92,19 @@ int main() {
universe.AddChild(std::move(outerMedium));
+ auto const allElementsInUniverse = std::invoke([&]() {
+ std::set<particles::Code> allElementsInUniverse;
+ auto collectElements = [&](auto& vtn) {
+ if (vtn.HasModelProperties()) {
+ auto const& comp =
+ vtn.GetModelProperties().GetNuclearComposition().GetComponents();
+ for (auto const c : comp) allElementsInUniverse.insert(c);
+ }
+ };
+ universe.walk(collectElements);
+ return allElementsInUniverse;
+ });
+
// setup particle stack, and add primary particle
setup::Stack stack;
stack.Clear();
@@ -132,7 +145,8 @@ int main() {
random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
random::RNGManager::GetInstance().RegisterRandomStream("pythia");
- process::sibyll::Interaction sibyll;
+
+ process::sibyll::Interaction sibyll(allElementsInUniverse);
process::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
process::sibyll::Decay decay;
process::particle_cut::ParticleCut cut(20_GeV);
diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc
index 7745c162..c8798021 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -85,6 +85,19 @@ int main() {
universe.AddChild(std::move(theMedium));
+ auto const allElementsInUniverse = std::invoke([&]() {
+ std::set<particles::Code> allElementsInUniverse;
+ auto collectElements = [&](auto& vtn) {
+ if (vtn.HasModelProperties()) {
+ auto const& comp =
+ vtn.GetModelProperties().GetNuclearComposition().GetComponents();
+ for (auto const c : comp) allElementsInUniverse.insert(c);
+ }
+ };
+ universe.walk(collectElements);
+ return allElementsInUniverse;
+ });
+
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
// setup particle stack, and add primary particle
@@ -126,7 +139,7 @@ int main() {
stack_inspector::StackInspector<setup::Stack> stackInspect(1, true, E0);
random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
- process::sibyll::Interaction sibyll;
+ process::sibyll::Interaction sibyll(allElementsInUniverse);
process::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
process::sibyll::Decay decay;
process::particle_cut::ParticleCut cut(20_GeV);
diff --git a/Processes/Sibyll/Interaction.cc b/Processes/Sibyll/Interaction.cc
index ec0bb142..6e7b9bb9 100644
--- a/Processes/Sibyll/Interaction.cc
+++ b/Processes/Sibyll/Interaction.cc
@@ -34,8 +34,6 @@ using Track = Trajectory;
namespace corsika::process::sibyll {
- Interaction::Interaction() {}
-
Interaction::~Interaction() {
cout << "Sibyll::Interaction n=" << fCount << " Nnuc=" << fNucCount << endl;
}
@@ -52,34 +50,19 @@ namespace corsika::process::sibyll {
}
// initialize cross sections
- InitializeCrossSectionTable();
+ InitializeCrossSectionTable(fTargetComposition);
}
- void Interaction::InitializeCrossSectionTable() {
+ void Interaction::InitializeCrossSectionTable(
+ std::set<corsika::particles::Code> allElementsInUniverse) {
using namespace units::si;
+ std::cout << "sibyll::Interaction: Initializing hadron-medium cross sections.."
+ << std::endl;
std::cout << "min,max:" << fLog10MinEnergyCoM << "," << fLog10MaxEnergyCoM
<< std::endl;
- // target nuclei
- // auto& universe = *(fEnvironment.GetUniverse());
- // auto const allElementsInUniverse = std::invoke([&]() {
- // std::set<particles::Code> allElementsInUniverse;
- // auto collectElements = [&](auto& vtn) {
- // if (vtn.HasModelProperties()) {
- // auto const& comp =
- // vtn.GetModelProperties().GetNuclearComposition().GetComponents();
- // for (auto const c : comp) allElementsInUniverse.insert(c);
- // }
- // };
- // universe.walk(collectElements);
- // return allElementsInUniverse;
- // });
-
- std::array<corsika::particles::Code, 2> allElementsInUniverse{
- corsika::particles::Code::Nitrogen, corsika::particles::Code::Oxygen};
-
for (auto& pproj : fProjectileHadrons) {
auto pproj_xcode = process::sibyll::GetSibyllXSCode(pproj);
for (auto& ptarg : allElementsInUniverse) {
@@ -102,6 +85,9 @@ namespace corsika::process::sibyll {
fHadronNuclearCrossSectionTables.insert(std::make_pair(col_conf, table));
}
}
+ std::cout << "sibyll::Interaction: cross sections for "
+ << fHadronNuclearCrossSectionTables.size() << " components initialized!"
+ << std::endl;
}
void Interaction::SetStable(std::vector<particles::Code> const& vParticleList) {
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
index 44e20286..13a02dcb 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/Processes/Sibyll/Interaction.h
@@ -15,6 +15,7 @@
#include <corsika/process/InteractionProcess.h>
#include <corsika/random/RNGManager.h>
#include <corsika/units/PhysicalUnits.h>
+#include <set>
#include <tuple>
namespace corsika::process::sibyll {
@@ -24,9 +25,15 @@ namespace corsika::process::sibyll {
int fCount = 0;
int fNucCount = 0;
bool fInitialized = false;
+ std::set<corsika::particles::Code> const fTargetComposition;
public:
- Interaction();
+ // Interaction();
+ // Interaction::Interaction()
+ // : fTargetComposition{corsika::particles::Code::Oxygen,
+ // corsika::particles::Code::Nitrogen} {}
+ Interaction(std::set<corsika::particles::Code> vComp)
+ : fTargetComposition(vComp) {}
~Interaction();
void Init();
@@ -56,7 +63,7 @@ namespace corsika::process::sibyll {
Energy range is 10GeV to 10PeV. Number of components in medium is arbitrary but
nucleon number is limited to A<18.
*/
- void InitializeCrossSectionTable();
+ void InitializeCrossSectionTable(std::set<corsika::particles::Code>);
/**
Calculate hadron nucleus/proton cross section in SIBYLL.
diff --git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc
index fa004958..0ab7187b 100644
--- a/Processes/Sibyll/testSibyll.cc
+++ b/Processes/Sibyll/testSibyll.cc
@@ -85,12 +85,13 @@ TEST_CASE("SibyllInterface", "[processes]") {
geometry::Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
1_km * std::numeric_limits<double>::infinity());
+ std::vector<particles::Code> elementsInMedium{particles::Code::Oxygen,
+ particles::Code::Helium};
+
using MyHomogeneousModel = environment::HomogeneousMedium<environment::IMediumModel>;
theMedium->SetModelProperties<MyHomogeneousModel>(
1_kg / (1_m * 1_m * 1_m),
- environment::NuclearComposition(
- std::vector<particles::Code>{particles::Code::Oxygen, particles::Code::Helium},
- std::vector<float>{0.5, 0.5}));
+ environment::NuclearComposition(elementsInMedium, std::vector<float>{0.5, 0.5}));
auto const* nodePtr = theMedium.get();
universe.AddChild(std::move(theMedium));
@@ -115,7 +116,8 @@ TEST_CASE("SibyllInterface", "[processes]") {
corsika::stack::SecondaryView view(particle);
auto projectile = view.GetProjectile();
- Interaction model;
+ Interaction model(
+ std::set<particles::Code>(elementsInMedium.begin(), elementsInMedium.end()));
model.Init();
[[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(projectile);
@@ -140,7 +142,8 @@ TEST_CASE("SibyllInterface", "[processes]") {
corsika::stack::SecondaryView view(particle);
auto projectile = view.GetProjectile();
- Interaction hmodel;
+ Interaction hmodel(
+ std::set<particles::Code>(elementsInMedium.begin(), elementsInMedium.end()));
NuclearInteraction model(hmodel, env);
model.Init();
--
GitLab