From 700f188cafdcd31e8ea9b413fac62cc0ff3e642e Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@kit.edu>
Date: Tue, 11 Dec 2018 20:12:06 +0100
Subject: [PATCH] homogeneous environment
---
Documentation/Examples/cascade_example.cc | 29 ++++++++++++++++-------
Environment/Environment.h | 10 ++++----
Environment/HomogeneousMedium.h | 6 ++---
Environment/NuclearComposition.h | 3 ++-
Environment/VolumeTreeNode.h | 8 +++----
Environment/testEnvironment.cc | 10 ++++++++
Framework/Cascade/testCascade.cc | 10 ++++----
Framework/Geometry/CoordinateSystem.cc | 3 ++-
Framework/Geometry/CoordinateSystem.h | 5 ++--
Processes/TrackingLine/TrackingLine.h | 26 +++++++++++++-------
10 files changed, 74 insertions(+), 36 deletions(-)
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index e4a8e5957..c3a8c9519 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -1,4 +1,3 @@
-
/**
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
@@ -18,6 +17,8 @@
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/environment/Environment.h>
+#include <corsika/environment/HomogeneousMedium.h>
+#include <corsika/environment/NuclearComposition.h>
#include <corsika/random/RNGManager.h>
@@ -29,6 +30,7 @@
#include <corsika/units/PhysicalUnits.h>
#include <iostream>
+#include <limits>
#include <typeinfo>
using namespace corsika;
@@ -192,9 +194,8 @@ public:
// rnd_ini_();
corsika::random::RNGManager& rmng = corsika::random::RNGManager::GetInstance();
- ;
- const std::string str_name = "s_rndm";
- rmng.RegisterRandomStream(str_name);
+
+ rmng.RegisterRandomStream("s_rndm");
// // corsika::random::RNG srng;
// auto srng = rmng.GetRandomStream("s_rndm");
@@ -242,13 +243,25 @@ double s_rndm_(int&) {
return rmng() / (double)rmng.max();
}
+class A {};
+class B : public A {};
+
int main() {
- corsika::environment::Environment env; // dummy environment
+ corsika::environment::Environment env; // dummy environment
auto& universe = *(env.GetUniverse());
auto theMedium = corsika::environment::Environment::CreateNode<Sphere>(
- Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m}, 100_km);
-
+ Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
+ 1_km * std::numeric_limits<double>::infinity());
+
+ using MyHomogeneousModel =
+ corsika::environment::HomogeneousMedium<corsika::environment::IMediumModel>;
+ theMedium->SetModelProperties<MyHomogeneousModel>(
+ 1_g / (1_m * 1_m * 1_m),
+ corsika::environment::NuclearComposition(
+ std::vector<corsika::particles::Code>{corsika::particles::Code::Proton},
+ std::vector<float>{1.}));
+
universe.AddChild(std::move(theMedium));
tracking_line::TrackingLine<setup::Stack> tracking(env);
@@ -261,7 +274,7 @@ int main() {
stack.Clear();
auto particle = stack.NewParticle();
- EnergyType E0 = 1_EeV;
+ EnergyType E0 = 100_TeV;
particle.SetEnergy(E0);
particle.SetPID(Code::Proton);
EAS.Init();
diff --git a/Environment/Environment.h b/Environment/Environment.h
index 285659a6b..5afe69c57 100644
--- a/Environment/Environment.h
+++ b/Environment/Environment.h
@@ -18,14 +18,16 @@
namespace corsika::environment {
struct Universe : public corsika::geometry::Volume {
- bool Contains(corsika::geometry::Point const&) const override {return true;}
+ bool Contains(corsika::geometry::Point const&) const override { return true; }
};
+ // template <typename IEnvironmentModel>
class Environment {
public:
- Environment() : fUniverse(std::make_unique<VolumeTreeNode<IEnvironmentModel>>(
- std::make_unique<Universe>())) {}
-
+ Environment()
+ : fUniverse(std::make_unique<VolumeTreeNode<IEnvironmentModel>>(
+ std::make_unique<Universe>())) {}
+
using IEnvironmentModel = corsika::setup::IEnvironmentModel;
auto& GetUniverse() { return fUniverse; }
diff --git a/Environment/HomogeneousMedium.h b/Environment/HomogeneousMedium.h
index bea6cd618..a8fa5f728 100644
--- a/Environment/HomogeneousMedium.h
+++ b/Environment/HomogeneousMedium.h
@@ -25,7 +25,7 @@
namespace corsika::environment {
template <class T>
- class HomogeneousMedium : T {
+ class HomogeneousMedium : public T {
corsika::units::si::MassDensityType const fDensity;
NuclearComposition const fNuclComp;
@@ -35,8 +35,8 @@ namespace corsika::environment {
: fDensity(pDensity)
, fNuclComp(pNuclComp){};
- corsika::units::si::MassDensityType GetMassDensity([
- [maybe_unused]] corsika::geometry::Point const& p) const override {
+ corsika::units::si::MassDensityType GetMassDensity(
+ corsika::geometry::Point const&) const override {
return fDensity;
}
NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
diff --git a/Environment/NuclearComposition.h b/Environment/NuclearComposition.h
index 484f98375..46a2bf5e8 100644
--- a/Environment/NuclearComposition.h
+++ b/Environment/NuclearComposition.h
@@ -3,6 +3,7 @@
#include <corsika/particles/ParticleProperties.h>
#include <numeric>
+#include <stdexcept>
#include <vector>
namespace corsika::environment {
@@ -20,7 +21,7 @@ namespace corsika::environment {
std::accumulate(pFractions.cbegin(), pFractions.cend(), 0.f);
if (!(0.999f < sumFractions && sumFractions < 1.001f)) {
- throw std::string("element fractions do not add up to 1");
+ throw std::runtime_error("element fractions do not add up to 1");
}
}
diff --git a/Environment/VolumeTreeNode.h b/Environment/VolumeTreeNode.h
index f5779f9ad..d92d77b39 100644
--- a/Environment/VolumeTreeNode.h
+++ b/Environment/VolumeTreeNode.h
@@ -87,12 +87,12 @@ namespace corsika::environment {
auto const& GetModelProperties() const { return *fModelProperties; }
- template <typename ModelProperties, typename... Args>
+ template <typename TModelProperties, typename... Args>
auto SetModelProperties(Args&&... args) {
- static_assert(std::is_base_of_v<IModelProperties, ModelProperties>,
+ static_assert(std::is_base_of_v<IModelProperties, TModelProperties>,
"unusable type provided");
- fModelProperties = std::make_shared<ModelProperties>(std::forward<Args>(args)...);
+ fModelProperties = std::make_shared<TModelProperties>(std::forward<Args>(args)...);
return fModelProperties;
}
@@ -111,7 +111,7 @@ namespace corsika::environment {
std::vector<VolumeTreeNode<IModelProperties> const*> fExcludedNodes;
VolumeTreeNode<IModelProperties> const* fParentNode = nullptr;
VolUPtr fGeoVolume;
- std::shared_ptr<IModelProperties> fModelProperties;
+ IMPSharedPtr fModelProperties;
};
} // namespace corsika::environment
diff --git a/Environment/testEnvironment.cc b/Environment/testEnvironment.cc
index 1efbc4e8f..8a8e9c13b 100644
--- a/Environment/testEnvironment.cc
+++ b/Environment/testEnvironment.cc
@@ -1,3 +1,13 @@
+/**
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
// cpp file
diff --git a/Framework/Cascade/testCascade.cc b/Framework/Cascade/testCascade.cc
index 1a72dc9ae..6bb0a10dd 100644
--- a/Framework/Cascade/testCascade.cc
+++ b/Framework/Cascade/testCascade.cc
@@ -1,4 +1,3 @@
-
/**
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
@@ -9,6 +8,8 @@
* the license.
*/
+#include <limits>
+
#include <corsika/environment/Environment.h>
#include <corsika/cascade/Cascade.h>
@@ -80,12 +81,13 @@ private:
TEST_CASE("Cascade", "[Cascade]") {
corsika::environment::Environment env; // dummy environment
auto& universe = *(env.GetUniverse());
- auto const radius = 20_km;
+ auto const radius = 1_m * std::numeric_limits<double>::infinity();
+ ;
auto theMedium = corsika::environment::Environment::CreateNode<Sphere>(
- Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m}, radius);
+ Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m}, radius);
universe.AddChild(std::move(theMedium));
- tracking_line::TrackingLine<setup::Stack> tracking(env);
+ tracking_line::TrackingLine<setup::Stack> tracking(env);
stack_inspector::StackInspector<setup::Stack> p0(true);
ProcessSplit p1;
diff --git a/Framework/Geometry/CoordinateSystem.cc b/Framework/Geometry/CoordinateSystem.cc
index 9cf3a1a64..cdc465c60 100644
--- a/Framework/Geometry/CoordinateSystem.cc
+++ b/Framework/Geometry/CoordinateSystem.cc
@@ -10,6 +10,7 @@
*/
#include <corsika/geometry/CoordinateSystem.h>
+#include <stdexcept>
using namespace corsika::geometry;
@@ -40,7 +41,7 @@ EigenTransform CoordinateSystem::GetTransformation(CoordinateSystem const& pFrom
commonBase = a;
} else {
- throw std::string("no connection between coordinate systems found!");
+ throw std::runtime_error("no connection between coordinate systems found!");
}
EigenTransform t = EigenTransform::Identity();
diff --git a/Framework/Geometry/CoordinateSystem.h b/Framework/Geometry/CoordinateSystem.h
index 3a06dfbd8..6f0584f5e 100644
--- a/Framework/Geometry/CoordinateSystem.h
+++ b/Framework/Geometry/CoordinateSystem.h
@@ -15,6 +15,7 @@
#include <corsika/geometry/QuantityVector.h>
#include <corsika/units/PhysicalUnits.h>
#include <Eigen/Dense>
+#include <stdexcept>
typedef Eigen::Transform<double, 3, Eigen::Affine> EigenTransform;
typedef Eigen::Translation<double, 3> EigenTranslation;
@@ -60,7 +61,7 @@ namespace corsika::geometry {
auto rotate(QuantityVector<phys::units::length_d> axis, double angle) const {
if (axis.eVector.isZero()) {
- throw std::string("null-vector given as axis parameter");
+ throw std::runtime_error("null-vector given as axis parameter");
}
EigenTransform const rotation{Eigen::AngleAxisd(angle, axis.eVector.normalized())};
@@ -71,7 +72,7 @@ namespace corsika::geometry {
auto translateAndRotate(QuantityVector<phys::units::length_d> translation,
QuantityVector<phys::units::length_d> axis, double angle) {
if (axis.eVector.isZero()) {
- throw std::string("null-vector given as axis parameter");
+ throw std::runtime_error("null-vector given as axis parameter");
}
EigenTransform const transf{Eigen::AngleAxisd(angle, axis.eVector.normalized()) *
diff --git a/Processes/TrackingLine/TrackingLine.h b/Processes/TrackingLine/TrackingLine.h
index 02c5e423c..b6b886408 100644
--- a/Processes/TrackingLine/TrackingLine.h
+++ b/Processes/TrackingLine/TrackingLine.h
@@ -24,6 +24,7 @@
#include <algorithm>
#include <iostream>
#include <optional>
+#include <stdexcept>
#include <utility>
using namespace corsika;
@@ -39,8 +40,9 @@ namespace corsika::process {
corsika::environment::Environment const& fEnvironment;
public:
- std::optional<std::pair<corsika::units::si::TimeType, corsika::units::si::TimeType>> TimeOfIntersection(
- corsika::geometry::Line const& traj, geometry::Sphere const& sphere) {
+ std::optional<std::pair<corsika::units::si::TimeType, corsika::units::si::TimeType>>
+ TimeOfIntersection(corsika::geometry::Line const& traj,
+ geometry::Sphere const& sphere) {
using namespace corsika::units::si;
auto const& cs = fEnvironment.GetCoordinateSystem();
geometry::Point const origin(cs, 0_m, 0_m, 0_m);
@@ -61,7 +63,8 @@ namespace corsika::process {
if (discriminant.magnitude() > 0) {
(-alpha - sqrt(discriminant)) / (2 * v0.squaredNorm());
- return std::make_pair((-alpha - sqrt(discriminant)) / (2 * v0.squaredNorm()), (-alpha + sqrt(discriminant)) / (2 * v0.squaredNorm()));
+ return std::make_pair((-alpha - sqrt(discriminant)) / (2 * v0.squaredNorm()),
+ (-alpha + sqrt(discriminant)) / (2 * v0.squaredNorm()));
} else {
return {};
}
@@ -92,13 +95,12 @@ namespace corsika::process {
volume); // for the moment we are a bit bold here and assume
// everything is a sphere, crashes with exception if not
- if (auto opt = TimeOfIntersection(line, sphere); opt.has_value())
- {
+ if (auto opt = TimeOfIntersection(line, sphere); opt.has_value()) {
auto const [t1, t2] = *opt;
if (t1.magnitude() >= 0)
- intersectionTimes.push_back(t1);
+ intersectionTimes.push_back(t1);
else if (t2.magnitude() >= 0)
- intersectionTimes.push_back(t2);
+ intersectionTimes.push_back(t2);
}
};
@@ -111,11 +113,17 @@ namespace corsika::process {
auto const minIter =
std::min_element(intersectionTimes.cbegin(), intersectionTimes.cend());
+ TimeType min;
+
if (minIter == intersectionTimes.cend()) {
- throw std::string("no intersection with anything!");
+ min = 1_s; // todo: do sth. more reasonable as soon as tracking is able
+ // to handle the numerics properly
+ //~ throw std::runtime_error("no intersection with anything!");
+ } else {
+ min = *minIter;
}
- return geometry::Trajectory<corsika::geometry::Line>(line, *minIter);
+ return geometry::Trajectory<corsika::geometry::Line>(line, min);
}
};
--
GitLab