From 90ff5ba8c8db6025e7094b75df00cf7468b6e526 Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@kit.edu>
Date: Mon, 14 Jan 2019 13:31:01 +0100
Subject: [PATCH] Cleanup
---
Documentation/Examples/cascade_example.cc | 2 +
Environment/Environment.h | 16 ++++---
Environment/HomogeneousMedium.h | 26 +-----------
Environment/IMediumModel.h | 23 +++++++++-
Environment/NuclearComposition.h | 3 ++
Framework/Cascade/Cascade.h | 1 +
Processes/Sibyll/Decay.h | 22 +++-------
Processes/Sibyll/Interaction.h | 52 ++++++++++++-----------
Processes/Sibyll/testSibyll.cc | 4 ++
9 files changed, 75 insertions(+), 74 deletions(-)
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 6dda88395..2ca239e5b 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -229,6 +229,8 @@ int main() {
// setup processes, decays and interactions
tracking_line::TrackingLine<setup::Stack> tracking(env);
stack_inspector::StackInspector<setup::Stack> p0(true);
+
+ corsika::random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
corsika::process::sibyll::Interaction sibyll(env);
corsika::process::sibyll::Decay decay;
ProcessCut cut(8_GeV);
diff --git a/Environment/Environment.h b/Environment/Environment.h
index aa567e876..1ab22e038 100644
--- a/Environment/Environment.h
+++ b/Environment/Environment.h
@@ -20,6 +20,8 @@
#include <limits>
namespace corsika::environment {
+ using BaseNodeType = VolumeTreeNode<corsika::setup::IEnvironmentModel>;
+
struct Universe : public corsika::geometry::Sphere {
Universe(corsika::geometry::CoordinateSystem const& pCS)
: corsika::geometry::Sphere(
@@ -37,7 +39,7 @@ namespace corsika::environment {
Environment()
: fCoordinateSystem{corsika::geometry::RootCoordinateSystem::GetInstance()
.GetRootCoordinateSystem()}
- , fUniverse(std::make_unique<VolumeTreeNode<IEnvironmentModel>>(
+ , fUniverse(std::make_unique<BaseNodeType>(
std::make_unique<Universe>(fCoordinateSystem))) {}
using IEnvironmentModel = corsika::setup::IEnvironmentModel;
@@ -48,19 +50,19 @@ namespace corsika::environment {
auto const& GetCoordinateSystem() const { return fCoordinateSystem; }
// factory method for creation of VolumeTreeNodes
- template <class VolumeType, typename... VolumeArgs>
- static auto CreateNode(VolumeArgs&&... args) {
- static_assert(std::is_base_of_v<corsika::geometry::Volume, VolumeType>,
+ template <class TVolumeType, typename... TVolumeArgs>
+ static auto CreateNode(TVolumeArgs&&... args) {
+ static_assert(std::is_base_of_v<corsika::geometry::Volume, TVolumeType>,
"unusable type provided, needs to be derived from "
"\"corsika::geometry::Volume\"");
- return std::make_unique<VolumeTreeNode<IEnvironmentModel>>(
- std::make_unique<VolumeType>(std::forward<VolumeArgs>(args)...));
+ return std::make_unique<BaseNodeType>(
+ std::make_unique<TVolumeType>(std::forward<TVolumeArgs>(args)...));
}
private:
corsika::geometry::CoordinateSystem const& fCoordinateSystem;
- VolumeTreeNode<IEnvironmentModel>::VTNUPtr fUniverse;
+ BaseNodeType::VTNUPtr fUniverse;
};
} // namespace corsika::environment
diff --git a/Environment/HomogeneousMedium.h b/Environment/HomogeneousMedium.h
index 9344d2009..d058f67c8 100644
--- a/Environment/HomogeneousMedium.h
+++ b/Environment/HomogeneousMedium.h
@@ -19,6 +19,8 @@
#include <corsika/random/RNGManager.h>
#include <corsika/units/PhysicalUnits.h>
+#include <cassert>
+
/**
* a homogeneous medium
*/
@@ -42,30 +44,6 @@ namespace corsika::environment {
}
NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
- corsika::particles::Code const& GetTarget(
- std::vector<corsika::units::si::CrossSectionType>& sigma) const override {
- using namespace corsika::units::si;
- int i = -1;
- corsika::units::si::CrossSectionType total_weighted_sigma = 0._mbarn;
- std::vector<float> fractions;
- for (auto w : fNuclComp.GetFractions()) {
- i++;
- std::cout << "HomogeneousMedium: fraction: " << w << std::endl;
- total_weighted_sigma += w * sigma[i];
- fractions.push_back(w * sigma[i] / 1_mbarn);
- }
-
- for (auto f : fractions) {
- f = f / (total_weighted_sigma / 1_mbarn);
- std::cout << "HomogeneousMedium: reweighted fraction: " << f << std::endl;
- }
- std::discrete_distribution channelDist(fractions.begin(), fractions.end());
- static corsika::random::RNG& kRNG =
- corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
- const int ichannel = channelDist(kRNG);
- return fNuclComp.GetComponents()[ichannel];
- }
-
corsika::units::si::GrammageType IntegratedGrammage(
corsika::geometry::Trajectory<corsika::geometry::Line> const&,
corsika::units::si::LengthType pTo) const override {
diff --git a/Environment/IMediumModel.h b/Environment/IMediumModel.h
index 3dabc3f98..4278a79a3 100644
--- a/Environment/IMediumModel.h
+++ b/Environment/IMediumModel.h
@@ -16,6 +16,7 @@
#include <corsika/geometry/Line.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/Trajectory.h>
+#include <random>
#include <corsika/units/PhysicalUnits.h>
namespace corsika::environment {
@@ -38,9 +39,27 @@ namespace corsika::environment {
corsika::units::si::GrammageType) const = 0;
virtual NuclearComposition const& GetNuclearComposition() const = 0;
+
+ template <class TRNG>
+ corsika::particles::Code SampleTarget(
+ std::vector<corsika::units::si::CrossSectionType> const& sigma, TRNG& randomStream) const {
+ using namespace corsika::units::si;
+
+ auto const& nuclComp = GetNuclearComposition();
+ auto const& fractions = nuclComp.GetFractions();
+ assert(sigma.size() == fractions.size());
+
+ std::vector<float> weights(fractions.size());
+
+ for (size_t i = 0; i < fractions.size(); ++i) {
+ std::cout << "HomogeneousMedium: fraction: " << fractions[i] << std::endl;
+ weights[i] = fractions[i] * sigma[i].magnitude();
+ }
- virtual corsika::particles::Code const& GetTarget(
- std::vector<corsika::units::si::CrossSectionType>&) const = 0;
+ std::discrete_distribution channelDist(weights.cbegin(), weights.cend());
+ const int iChannel = channelDist(randomStream);
+ return nuclComp.GetComponents()[iChannel];
+ }
};
} // namespace corsika::environment
diff --git a/Environment/NuclearComposition.h b/Environment/NuclearComposition.h
index a71f22305..3bc1101a4 100644
--- a/Environment/NuclearComposition.h
+++ b/Environment/NuclearComposition.h
@@ -15,6 +15,7 @@
#include <corsika/particles/ParticleProperties.h>
#include <numeric>
#include <stdexcept>
+#include <cassert>
#include <vector>
namespace corsika::environment {
@@ -28,6 +29,8 @@ namespace corsika::environment {
std::vector<float> pFractions)
: fNumberFractions(pFractions)
, fComponents(pComponents) {
+
+ assert(pComponents.size() == pFractions.size());
auto const sumFractions =
std::accumulate(pFractions.cbegin(), pFractions.cend(), 0.f);
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index ba05bd78e..999894f70 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -56,6 +56,7 @@ namespace corsika::cascade {
}
}
+private:
void Step(Particle& particle) {
using namespace corsika::units::si;
diff --git a/Processes/Sibyll/Decay.h b/Processes/Sibyll/Decay.h
index e7ed3165d..217ed87fa 100644
--- a/Processes/Sibyll/Decay.h
+++ b/Processes/Sibyll/Decay.h
@@ -21,8 +21,6 @@
#include <corsika/particles/ParticleProperties.h>
-#include <fenv.h>
-
namespace corsika::process {
namespace sibyll {
@@ -81,7 +79,7 @@ namespace corsika::process {
// i.e. corsika::particles::ListOfParticles()
for (auto& p : corsika2sibyll) {
// std::cout << (int)p << std::endl;
- const int sibCode = (int)p;
+ const int sibCode = static_cast<int>(p);
// skip unknown and antiparticles
if (sibCode < 1) continue;
s_csydec_.idb[sibCode - 1] = -1 * abs(s_csydec_.idb[sibCode - 1]);
@@ -114,7 +112,7 @@ namespace corsika::process {
}
// set Leptons and Proton and Neutron stable
// use stack to loop over particles
- const std::vector<corsika::particles::Code> particleList = {
+ constexpr corsika::particles::Code particleList[] = {
corsika::particles::Code::Proton, corsika::particles::Code::Neutron,
corsika::particles::Code::Electron, corsika::particles::Code::Positron,
corsika::particles::Code::NuE, corsika::particles::Code::NuEBar,
@@ -131,7 +129,7 @@ namespace corsika::process {
}
template <typename Particle>
- corsika::units::si::TimeType GetLifetime(Particle& p) {
+ corsika::units::si::TimeType GetLifetime(Particle const& p) {
using std::cout;
using std::endl;
using namespace corsika::units::si;
@@ -169,11 +167,6 @@ namespace corsika::process {
using corsika::geometry::Point;
using namespace corsika::units::si;
- // TODO: this should be done in a central, common place. Not here..
-#ifndef CORSIKA_OSX
- feenableexcept(FE_INVALID);
-#endif
-
fCount++;
SibStack ss;
ss.Clear();
@@ -189,8 +182,8 @@ namespace corsika::process {
// with sibyll internal values
pin.SetMass(corsika::particles::GetMass(pCode));
// remember position
- Point decayPoint = p.GetPosition();
- TimeType t0 = p.GetTime();
+ Point const decayPoint = p.GetPosition();
+ TimeType const t0 = p.GetTime();
// remove original particle from corsika stack
p.Delete();
// set all particles/hadrons unstable
@@ -219,11 +212,6 @@ namespace corsika::process {
}
// empty sibyll stack
ss.Clear();
-
- // TODO: this should be done in a central, common place. Not here..
-#ifndef CORSIKA_OSX
- fedisableexcept(FE_INVALID);
-#endif
}
};
} // namespace sibyll
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
index 003d22a5f..f35b244f8 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/Processes/Sibyll/Interaction.h
@@ -45,16 +45,13 @@ namespace corsika::process::sibyll {
using std::cout;
using std::endl;
- corsika::random::RNGManager& rmng = corsika::random::RNGManager::GetInstance();
- rmng.RegisterRandomStream("s_rndm");
-
// initialize Sibyll
sibyll_ini_();
}
std::tuple<corsika::units::si::CrossSectionType, int> GetCrossSection(
- const corsika::particles::Code& BeamId, const corsika::particles::Code& TargetId,
- const corsika::units::si::HEPEnergyType& CoMenergy) {
+ const corsika::particles::Code BeamId, const corsika::particles::Code TargetId,
+ const corsika::units::si::HEPEnergyType CoMenergy) {
using namespace corsika::units::si;
double sigProd, dummy, dum1, dum2, dum3, dum4;
double dumdif[3];
@@ -67,19 +64,18 @@ namespace corsika::process::sibyll {
"Sibyll target outside range. Only nuclei with A<18 are allowed.");
sib_sigma_hnuc_(iBeam, iTarget, dEcm, sigProd, dummy);
return std::make_tuple(sigProd * 1_mbarn, iTarget);
+ } else if (TargetId == corsika::particles::Proton::GetCode()) {
+ sib_sigma_hp_(iBeam, dEcm, dum1, dum2, sigProd, dumdif, dum3, dum4);
+ return std::make_tuple(sigProd * 1_mbarn, 1);
} else {
- if (TargetId == corsika::particles::Proton::GetCode()) {
- sib_sigma_hp_(iBeam, dEcm, dum1, dum2, sigProd, dumdif, dum3, dum4);
- return std::make_tuple(sigProd * 1_mbarn, 1);
- } else
- // no interaction in sibyll possible, return infinite cross section? or throw?
- sigProd = std::numeric_limits<double>::infinity();
+ // no interaction in sibyll possible, return infinite cross section? or throw?
+ sigProd = std::numeric_limits<double>::infinity();
return std::make_tuple(sigProd * 1_mbarn, 0);
}
}
template <typename Particle, typename Track>
- corsika::units::si::GrammageType GetInteractionLength(Particle& p, Track&) {
+ corsika::units::si::GrammageType GetInteractionLength(Particle const& p, Track&) {
using namespace corsika::units;
using namespace corsika::units::si;
@@ -105,11 +101,13 @@ namespace corsika::process::sibyll {
// total momentum and energy
HEPEnergyType Elab = p.GetEnergy() + nucleon_mass;
- MomentumVector Ptot(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
- Ptot += p.GetMomentum();
- Ptot += pTarget;
+ MomentumVector pTotLab(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
+ pTotLab += p.GetMomentum();
+ pTotLab += pTarget;
+ auto const pTotLabNorm = pTotLab.norm();
// calculate cm. energy
- const HEPEnergyType ECoM = sqrt(Elab * Elab - Ptot.squaredNorm());
+ const HEPEnergyType ECoM = sqrt(
+ (Elab + pTotLabNorm) * (Elab - pTotLabNorm)); // binomial for numerical accuracy
std::cout << "Interaction: LambdaInt: \n"
<< " input energy: " << p.GetEnergy() / 1_GeV << endl
@@ -136,7 +134,7 @@ namespace corsika::process::sibyll {
// get weights of components from environment/medium
const auto w = mediumComposition.GetFractions();
// loop over components in medium
- for (auto targetId : mediumComposition.GetComponents()) {
+ for (auto const targetId : mediumComposition.GetComponents()) {
i++;
cout << "Interaction: get interaction length for target: " << targetId << endl;
@@ -155,7 +153,7 @@ namespace corsika::process::sibyll {
// calculate interaction length in medium
#warning check interaction length units
- GrammageType int_length =
+ GrammageType const int_length =
avgTargetMassNumber * corsika::units::constants::u / weightedProdCrossSection;
std::cout << "Interaction: "
<< "interaction length (g/cm2): "
@@ -238,7 +236,7 @@ namespace corsika::process::sibyll {
// sample target mass number
const auto currentNode = fEnvironment.GetUniverse()->GetContainingNode(pOrig);
- const auto mediumComposition =
+ const auto& mediumComposition =
currentNode->GetModelProperties().GetNuclearComposition();
// get cross sections for target materials
/*
@@ -246,14 +244,18 @@ namespace corsika::process::sibyll {
should be passed from GetInteractionLength if possible
*/
#warning reading interaction cross section again, should not be necessary
- std::vector<si::CrossSectionType> cross_section_of_components;
- for (auto targetId : mediumComposition.GetComponents()) {
+ auto const& compVec = mediumComposition.GetComponents();
+ std::vector<si::CrossSectionType> cross_section_of_components(
+ mediumComposition.GetComponents().size());
+
+ for (size_t i = 0; i < compVec.size(); ++i) {
+ auto const targetId = compVec[i];
const auto [sigProd, nNuc] = GetCrossSection(corsikaBeamId, targetId, Ecm);
- cross_section_of_components.push_back(sigProd);
+ cross_section_of_components[i] = sigProd;
}
- const auto targetCode =
- currentNode->GetModelProperties().GetTarget(cross_section_of_components);
+ const auto targetCode = currentNode->GetModelProperties().SampleTarget(
+ cross_section_of_components, fRNG);
cout << "Interaction: target selected: " << targetCode << endl;
/*
FOR NOW: allow nuclei with A<18 or protons only.
@@ -336,6 +338,8 @@ namespace corsika::process::sibyll {
private:
corsika::environment::Environment const& fEnvironment;
+ corsika::random::RNG& fRNG =
+ corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
};
} // namespace corsika::process::sibyll
diff --git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc
index 8c8ce3b5e..9b078fd08 100644
--- a/Processes/Sibyll/testSibyll.cc
+++ b/Processes/Sibyll/testSibyll.cc
@@ -13,6 +13,8 @@
#include <corsika/process/sibyll/Interaction.h>
#include <corsika/process/sibyll/ParticleConversion.h>
+#include <corsika/random/RNGManager.h>
+
#include <corsika/particles/ParticleProperties.h>
#include <corsika/geometry/Point.h>
@@ -109,6 +111,8 @@ TEST_CASE("SibyllInterface", "[processes]") {
geometry::Line line(origin, v);
geometry::Trajectory<geometry::Line> track(line, 10_s);
+ corsika::random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
+
SECTION("InteractionInterface") {
setup::Stack stack;
--
GitLab