IAP GITLAB

Skip to content
Snippets Groups Projects
Commit 0420d736 authored by ralfulrich's avatar ralfulrich
Browse files

a few merge conflicts

parents beca7928 c2686551
No related branches found
No related tags found
1 merge request!27Resolve "Process selection logic"
Pipeline #28 failed
Showing
with 512 additions and 32 deletions
......@@ -65,6 +65,7 @@ find_package (Eigen3 REQUIRED)
add_subdirectory (ThirdParty)
#add_subdirectory (Utilities)
add_subdirectory (Framework)
add_subdirectory (Environment)
add_subdirectory (Stack)
add_subdirectory (Setup)
add_subdirectory (Processes)
......
/**
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
......@@ -17,22 +16,33 @@
#include <corsika/setup/SetupStack.h>
#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>
#include <corsika/cascade/SibStack.h>
#include <corsika/cascade/sibyll2.3c.h>
#include <corsika/geometry/Sphere.h>
#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/units/PhysicalUnits.h>
#include <iostream>
#include <limits>
#include <typeinfo>
using namespace corsika;
using namespace corsika::process;
using namespace corsika::units;
using namespace corsika::particles;
using namespace corsika::random;
using namespace corsika::geometry;
using namespace corsika::environment;
#include <iostream>
#include <typeinfo>
using namespace std;
using namespace corsika::units::si;
static int fCount = 0;
......@@ -95,7 +105,6 @@ public:
template <typename Particle, typename Track, typename Stack>
EProcessReturn DoContinuous(Particle&, Track&, Stack&) const {
// corsika::utls::ignore(p);
return EProcessReturn::eOk;
}
......@@ -187,9 +196,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");
......@@ -237,9 +245,28 @@ double s_rndm_(int&) {
return rmng() / (double)rmng.max();
}
class A {};
class B : public A {};
int main() {
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},
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;
tracking_line::TrackingLine<setup::Stack> tracking(env);
stack_inspector::StackInspector<setup::Stack> p0(true);
ProcessSplit p1;
const auto sequence = p0 + p1;
......@@ -249,7 +276,7 @@ int main() {
stack.Clear();
auto particle = stack.NewParticle();
EnergyType E0 = 100_GeV;
EnergyType E0 = 100_TeV;
particle.SetEnergy(E0);
particle.SetPID(Code::Proton);
EAS.Init();
......
set (
ENVIRONMENT_HEADERS
VolumeTreeNode.h
IMediumModel.h
NuclearComposition.h
HomogeneousMedium.h
Environment.h
)
set (
ENVIRONMENT_NAMESPACE
corsika/environment
)
add_library (CORSIKAenvironment INTERFACE)
CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAenvironment ${ENVIRONMENT_NAMESPACE} ${ENVIRONMENT_HEADERS})
# target dependencies on other libraries (also the header onlys)
target_link_libraries (
CORSIKAenvironment
INTERFACE
CORSIKAgeometry
CORSIKAparticles
CORSIKAunits
)
target_include_directories (
CORSIKAenvironment
INTERFACE
$<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
$<INSTALL_INTERFACE:include>
)
install (
TARGETS CORSIKAenvironment
LIBRARY DESTINATION lib
ARCHIVE DESTINATION lib
PUBLIC_HEADER DESTINATION include/${ENVIRONMENT_NAMESPACE}
)
# --------------------
# code unit testing
add_executable (testEnvironment testEnvironment.cc)
target_link_libraries (
testEnvironment
CORSIKAenvironment
CORSIKAthirdparty # for catch2
)
add_test (NAME testGeometry COMMAND testGeometry)
/**
* (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.
*/
#ifndef _include_Environment_h
#define _include_Environment_h
#include <corsika/environment/IMediumModel.h>
#include <corsika/environment/VolumeTreeNode.h>
#include <corsika/geometry/RootCoordinateSystem.h>
#include <corsika/setup/SetupEnvironment.h>
namespace corsika::environment {
struct Universe : public corsika::geometry::Volume {
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>())) {}
using IEnvironmentModel = corsika::setup::IEnvironmentModel;
auto& GetUniverse() { return fUniverse; }
auto const& GetUniverse() const { return fUniverse; }
auto const& GetCoordinateSystem() const {
return corsika::geometry::RootCoordinateSystem::GetInstance().GetRootCS();
}
// factory method for creation of VolumeTreeNodes
template <class VolumeType, typename... Args>
static auto CreateNode(Args&&... args) {
static_assert(std::is_base_of_v<corsika::geometry::Volume, VolumeType>,
"unusable type provided, needs to be derived from "
"\"corsika::geometry::Volume\"");
return std::make_unique<VolumeTreeNode<IEnvironmentModel>>(
std::make_unique<VolumeType>(std::forward<Args>(args)...));
}
private:
VolumeTreeNode<IEnvironmentModel>::VTNUPtr fUniverse;
};
} // namespace corsika::environment
#endif
/**
* (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.
*/
#ifndef _include_HomogeneousMedium_h_
#define _include_HomogeneousMedium_h_
#include <corsika/environment/NuclearComposition.h>
#include <corsika/geometry/Line.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/Trajectory.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/units/PhysicalUnits.h>
/**
* a homogeneous medium
*/
namespace corsika::environment {
template <class T>
class HomogeneousMedium : public T {
corsika::units::si::MassDensityType const fDensity;
NuclearComposition const fNuclComp;
public:
HomogeneousMedium(corsika::units::si::MassDensityType pDensity,
NuclearComposition pNuclComp)
: fDensity(pDensity)
, fNuclComp(pNuclComp){};
corsika::units::si::MassDensityType GetMassDensity(
corsika::geometry::Point const&) const override {
return fDensity;
}
NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
corsika::units::si::GrammageType IntegratedGrammage(
corsika::geometry::Trajectory<corsika::geometry::Line> const& pTraj,
corsika::units::si::TimeType pTo) const override {
using namespace corsika::units::si;
return pTraj.ArcLength(0_s, pTo) * fDensity;
}
corsika::units::si::TimeType FromGrammage(
corsika::geometry::Trajectory<corsika::geometry::Line> const& pTraj,
corsika::units::si::GrammageType pGrammage) const override {
return pTraj.TimeFromArclength(pGrammage / fDensity);
}
};
} // namespace corsika::environment
#endif
#ifndef _include_IMediumModel_h
#define _include_IMediumModel_h
#include <corsika/environment/NuclearComposition.h>
#include <corsika/geometry/Line.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/Trajectory.h>
#include <corsika/units/PhysicalUnits.h>
namespace corsika::environment {
class IMediumModel {
public:
virtual ~IMediumModel() = default;
virtual corsika::units::si::MassDensityType GetMassDensity(
corsika::geometry::Point const&) const = 0;
// todo: think about the mixin inheritance of the trajectory vs the BaseTrajectory
// approach for now, only lines are supported
virtual corsika::units::si::GrammageType IntegratedGrammage(
corsika::geometry::Trajectory<corsika::geometry::Line> const&,
corsika::units::si::TimeType) const = 0;
virtual corsika::units::si::TimeType FromGrammage(
corsika::geometry::Trajectory<corsika::geometry::Line> const&,
corsika::units::si::GrammageType) const = 0;
virtual NuclearComposition const& GetNuclearComposition() const = 0;
};
} // namespace corsika::environment
#endif
#ifndef _include_NuclearComposition_h
#define _include_NuclearComposition_h
#include <corsika/particles/ParticleProperties.h>
#include <numeric>
#include <stdexcept>
#include <vector>
namespace corsika::environment {
class NuclearComposition {
std::vector<float> const fNumberFractions; //!< relative fractions of number density
std::vector<corsika::particles::Code> const
fComponents; //!< particle codes of consitutents
public:
NuclearComposition(std::vector<corsika::particles::Code> pComponents,
std::vector<float> pFractions)
: fNumberFractions(pFractions)
, fComponents(pComponents) {
auto const sumFractions =
std::accumulate(pFractions.cbegin(), pFractions.cend(), 0.f);
if (!(0.999f < sumFractions && sumFractions < 1.001f)) {
throw std::runtime_error("element fractions do not add up to 1");
}
}
auto size() const { return fNumberFractions.size(); }
auto const& GetFractions() const { return fNumberFractions; }
auto const& GetComponents() const { return fComponents; }
};
} // namespace corsika::environment
#endif
/**
* (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.
*/
#ifndef _include_VolumeTreeNode_H
#define _include_VolumeTreeNode_H
#include <corsika/geometry/Volume.h>
#include <memory>
#include <vector>
namespace corsika::environment {
class Empty {}; //<! intended for usage as default template argument
template <typename IModelProperties = Empty>
class VolumeTreeNode {
public:
using VTNUPtr = std::unique_ptr<VolumeTreeNode<IModelProperties>>;
using IMPSharedPtr = std::shared_ptr<IModelProperties>;
using VolUPtr = std::unique_ptr<corsika::geometry::Volume>;
VolumeTreeNode(VolUPtr pVolume = nullptr)
: fGeoVolume(std::move(pVolume)) {}
bool Contains(corsika::geometry::Point const& p) const {
return fGeoVolume->Contains(p);
}
VolumeTreeNode<IModelProperties> const* Excludes(
corsika::geometry::Point const& p) const {
auto exclContainsIter =
std::find_if(fExcludedNodes.cbegin(), fExcludedNodes.cend(),
[&](auto const& s) { return bool(s->Contains(p)); });
return exclContainsIter != fExcludedNodes.cend() ? *exclContainsIter : nullptr;
}
/** returns a pointer to the sub-VolumeTreeNode which is "responsible" for the given
* \class Point \arg p, or nullptr iff \arg p is not contained in this volume.
*/
VolumeTreeNode<IModelProperties> const* GetContainingNode(
corsika::geometry::Point const& p) const {
if (!Contains(p)) { return nullptr; }
if (auto const childContainsIter =
std::find_if(fChildNodes.cbegin(), fChildNodes.cend(),
[&](auto const& s) { return bool(s->Contains(p)); });
childContainsIter == fChildNodes.cend()) // not contained in any of the children
{
if (auto const exclContainsIter = Excludes(p)) // contained in any excluded nodes
{
return exclContainsIter->GetContainingNode(p);
} else {
return this;
}
} else {
return (*childContainsIter)->GetContainingNode(p);
}
}
void AddChild(VTNUPtr pChild) {
pChild->fParentNode = this;
fChildNodes.push_back(std::move(pChild));
// It is a bad idea to return an iterator to the inserted element
// because it might get invalidated when the vector needs to grow
// later and the caller won't notice.
}
void ExcludeOverlapWith(VTNUPtr const& pNode) {
fExcludedNodes.push_back(pNode.get());
}
auto* GetParent() const { return fParentNode; };
auto const& GetChildNodes() const { return fChildNodes; }
auto const& GetExcludedNodes() const { return fExcludedNodes; }
auto const& GetVolume() const { return *fGeoVolume; }
auto const& GetModelProperties() const { return *fModelProperties; }
template <typename TModelProperties, typename... Args>
auto SetModelProperties(Args&&... args) {
static_assert(std::is_base_of_v<IModelProperties, TModelProperties>,
"unusable type provided");
fModelProperties = std::make_shared<TModelProperties>(std::forward<Args>(args)...);
return fModelProperties;
}
void SetModelProperties(IMPSharedPtr ptr) { fModelProperties = ptr; }
template <class MediumType, typename... Args>
static auto CreateMedium(Args&&... args) {
static_assert(std::is_base_of_v<IMediumModel, MediumType>,
"unusable type provided, needs to be derived from \"IMediumModel\"");
return std::make_shared<MediumType>(std::forward<Args>(args)...);
}
private:
std::vector<VTNUPtr> fChildNodes;
std::vector<VolumeTreeNode<IModelProperties> const*> fExcludedNodes;
VolumeTreeNode<IModelProperties> const* fParentNode = nullptr;
VolUPtr fGeoVolume;
IMPSharedPtr fModelProperties;
};
} // namespace corsika::environment
#endif
/**
* (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
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/IMediumModel.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/environment/VolumeTreeNode.h>
#include <corsika/particles/ParticleProperties.h>
#include <catch2/catch.hpp>
using namespace corsika::geometry;
using namespace corsika::environment;
using namespace corsika::units::si;
TEST_CASE("HomogeneousMedium") {
NuclearComposition const protonComposition(
std::vector<corsika::particles::Code>{corsika::particles::Code::Proton},
std::vector<float>{1.f});
HomogeneousMedium<IMediumModel> const medium(19.2_g / cube(1_cm), protonComposition);
}
/**
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
......@@ -9,6 +8,10 @@
* the license.
*/
#include <limits>
#include <corsika/environment/Environment.h>
#include <corsika/cascade/Cascade.h>
#include <corsika/process/ProcessSequence.h>
......@@ -38,6 +41,7 @@ using namespace corsika::geometry;
#include <iostream>
using namespace std;
using namespace corsika::units::si;
static int fCount = 0;
......@@ -82,7 +86,15 @@ TEST_CASE("Cascade", "[Cascade]") {
const std::string str_name = "s_rndm";
rmng.RegisterRandomStream(str_name);
tracking_line::TrackingLine<setup::Stack> tracking;
corsika::environment::Environment env; // dummy environment
auto& universe = *(env.GetUniverse());
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);
universe.AddChild(std::move(theMedium));
tracking_line::TrackingLine<setup::Stack> tracking(env);
stack_inspector::StackInspector<setup::Stack> p0(true);
ProcessSplit p1;
......
......@@ -42,7 +42,10 @@ namespace corsika::geometry {
* parameterized by \arg t1 and \arg t2. Requires \arg t2 > \arg t1.
*/
virtual LengthType GetDistance(corsika::units::si::TimeType t1,
virtual corsika::units::si::TimeType TimeFromArclength(
corsika::units::si::LengthType) const = 0;
virtual LengthType ArcLength(corsika::units::si::TimeType t1,
corsika::units::si::TimeType t2) const = 0;
virtual corsika::units::si::TimeType GetDuration(
......
......@@ -10,6 +10,7 @@ set (
Point.h
Line.h
Sphere.h
Volume.h
CoordinateSystem.h
RootCoordinateSystem.h
Helix.h
......
......@@ -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();
......
......@@ -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()) *
......
......@@ -58,10 +58,15 @@ namespace corsika::geometry {
auto GetRadius() const { return radius; }
LengthType GetDistanceBetween(corsika::units::si::TimeType t1,
corsika::units::si::TimeType t2) const {
corsika::units::si::LengthType ArcLength(
corsika::units::si::TimeType t1, corsika::units::si::TimeType t2) const {
return (vPar + vPerp).norm() * (t2 - t1);
}
corsika::units::si::TimeType TimeFromArclength(
corsika::units::si::LengthType t) const {
return t / (vPar + vPerp).norm();
}
};
} // namespace corsika::geometry
......
......@@ -32,11 +32,18 @@ namespace corsika::geometry {
Point GetPosition(corsika::units::si::TimeType t) const { return r0 + v0 * t; }
LengthType GetDistanceBetween(corsika::units::si::TimeType t1,
corsika::units::si::TimeType t2) const {
// assert(t2 >= t1);
LengthType ArcLength(corsika::units::si::TimeType t1,
corsika::units::si::TimeType t2) const {
return v0.norm() * (t2 - t1);
}
corsika::units::si::TimeType TimeFromArclength(
corsika::units::si::LengthType t) const {
return t / v0.norm();
}
auto GetR0() const { return r0; }
auto GetV0() const { return v0; }
};
} // namespace corsika::geometry
......
......@@ -13,23 +13,27 @@
#define _include_SPHERE_H_
#include <corsika/geometry/Point.h>
#include <corsika/geometry/Volume.h>
#include <corsika/units/PhysicalUnits.h>
namespace corsika::geometry {
class Sphere {
Point center;
LengthType const radius;
class Sphere : public Volume {
Point const fCenter;
LengthType const fRadius;
public:
Sphere(Point const& pCenter, LengthType const pRadius)
: center(pCenter)
, radius(pRadius) {}
: fCenter(pCenter)
, fRadius(pRadius) {}
//! returns true if the Point \a p is within the sphere
auto Contains(Point const& p) const {
return radius * radius > (center - p).squaredNorm();
//! returns true if the Point p is within the sphere
bool Contains(Point const& p) const override {
return fRadius * fRadius > (fCenter - p).squaredNorm();
}
auto& GetCenter() const { return fCenter; }
auto GetRadius() const { return fRadius; }
};
} // namespace corsika::geometry
......
......@@ -13,6 +13,7 @@
#define _include_TRAJECTORY_H
#include <corsika/units/PhysicalUnits.h>
#include <corsika/geometry/Point.h>
using corsika::units::si::LengthType;
using corsika::units::si::TimeType;
......@@ -25,7 +26,7 @@ namespace corsika::geometry {
corsika::units::si::TimeType fTimeLength;
public:
using T::GetDistanceBetween;
using T::ArcLength;
using T::GetPosition;
Trajectory(T const& theT, corsika::units::si::TimeType timeLength)
......@@ -36,14 +37,14 @@ namespace corsika::geometry {
return fTraj.GetPosition(t + fTStart);
}*/
Point GetPosition(const double u) const { return T::GetPosition(fTimeLength * u); }
Point GetPosition(double u) const { return T::GetPosition(fTimeLength * u); }
TimeType GetDuration() const { return fTimeLength; }
LengthType GetDistance(const corsika::units::si::TimeType t) const {
LengthType GetDistance(corsika::units::si::TimeType t) const {
assert(t > fTimeLength);
assert(t >= 0 * corsika::units::si::second);
return T::DistanceBetween(0, t);
return T::ArcLength(0, t);
}
};
......
......@@ -181,6 +181,17 @@ namespace corsika::geometry {
bareResult);
}
}
template <typename dim2>
auto dot(Vector<dim2> pV) const {
auto const c1 = GetComponents().eVector;
auto const c2 = pV.GetComponents(*BaseVector<dim>::cs).eVector;
auto const bareResult = c1.dot(c2);
using ProdQuantity = phys::units::detail::Product<dim, dim2, double, double>;
return ProdQuantity(phys::units::detail::magnitude_tag, bareResult);
}
};
} // namespace corsika::geometry
......
#ifndef _include_VOLUME_H_
#define _include_VOLUME_H_
#include <corsika/geometry/Point.h>
namespace corsika::geometry {
class Volume {
public:
//! returns true if the Point p is within the volume
virtual bool Contains(Point const& p) const = 0;
virtual ~Volume() = default;
};
} // namespace corsika::geometry
#endif
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment