diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc index e4a8e5957044bbf2614f66a8c0dcc2aadff3cae5..c3a8c95191effe7928ac470281617c4f34d486dd 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 285659a6bc1ec707673f55ee03f1f8d9ab1a9ea1..5afe69c57000976ef529559696e6bcdca31da609 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 bea6cd6181066697636bacd4ca34abec5391bda2..a8fa5f7283cfe0e776383f24915b98748d54349f 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 484f9837547c13996da409beb70d99a589cfc688..46a2bf5e859f10a7f5f13a4e4312d38cb1ae6086 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 f5779f9ad339528e56ba4a5f14078af5afcca8e2..d92d77b398bd180d9d64dc5039f32e50ca9be487 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 1efbc4e8f62bd90e6b2a9a53a966f3da890cf32b..8a8e9c13b5a899ad983d0e852379e5621c9f1489 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 1a72dc9ae5e14559670280253f8f301e9bdec9b4..6bb0a10ddd39956df11e469739faf4dce3cfef7a 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 9cf3a1a64125d4fe53236529bd821c6b2c5e57c5..cdc465c60842fb99d2f73539dd657d63131ef251 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 3a06dfbd8bc65b06a83a58184088d5bec42e8b03..6f0584f5e015aa521a0031919e0b0dda98d8d224 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 02c5e423c2a7ef2b52d31253626086c78c5b0a3f..b6b8864087b53fe577029c379996ca9d6ef29d35 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); } };