diff --git a/CMakeLists.txt b/CMakeLists.txt
index 1ac8449f081cd503c6c276eff57591562f34c2fd..b6969bead512aacffc631f993a6a77d6e4f06bff 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -39,6 +39,10 @@ if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
"Debug" "Release" "MinSizeRel" "RelWithDebInfo")
endif()
+#set(CMAKE_CXX_FLAGS "-Wall -Wextra")
+set(CMAKE_CXX_FLAGS_DEBUG "-O0 -g")
+set(CMAKE_CXX_FLAGS_RELEASE "-O3")
+
# unit testing coverage, does not work yet
#include (CodeCoverage)
##set(COVERAGE_LCOV_EXCLUDES 'Documentation/*')
@@ -62,6 +66,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)
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index e6e2fe483227a826a1b677be1212dbb4648b1d12..f9f344d4470a9a5638a00d8c041c87eb90a0b27b 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
*
@@ -17,26 +16,36 @@
#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/process/sibyll/ProcessDecay.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::setup;
+using namespace corsika::geometry;
+using namespace corsika::environment;
-#include <iostream>
-#include <typeinfo>
using namespace std;
+using namespace corsika::units::si;
static int fCount = 0;
static EnergyType fEnergy = 0. * 1_GeV;
@@ -197,7 +206,7 @@ public:
private:
};
-class ProcessSplit : public corsika::process::BaseProcess<ProcessSplit> {
+class ProcessSplit : public corsika::process::InteractionProcess<ProcessSplit> {
public:
ProcessSplit() {}
@@ -226,8 +235,8 @@ public:
}
}
- template <typename Particle>
- double MinStepLength(Particle& p, setup::Trajectory&) const {
+ template <typename Particle, typename Track>
+ double GetInteractionLength(Particle& p, Track&) const {
// coordinate system, get global frame of reference
CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
@@ -239,8 +248,8 @@ public:
int kBeam = process::sibyll::GetSibyllXSCode(corsikaBeamId);
bool kInteraction = process::sibyll::CanInteract(corsikaBeamId);
-
- /*
+
+ /*
the target should be defined by the Environment,
ideally as full particle object so that the four momenta
and the boosts can be defined..
@@ -268,7 +277,7 @@ public:
<< " beam pid:" << p.GetPID() << endl
<< " target mass number:" << kTarget << std::endl;
- double next_step;
+ double int_length = 0;
if (kInteraction) {
double prodCrossSection, dummy, dum1, dum2, dum3, dum4;
@@ -289,35 +298,35 @@ public:
std::cout << "ProcessSplit: "
<< "nucleon mass " << nucleon_mass << std::endl;
// calculate interaction length in medium
- double int_length = kTarget * (nucleon_mass / 1_g) / (sig / 1_cmeter / 1_cmeter);
+ int_length = kTarget * (nucleon_mass / 1_g) / (sig / 1_cmeter / 1_cmeter);
// pick random step lenth
std::cout << "ProcessSplit: "
<< "interaction length (g/cm2): " << int_length << std::endl;
- // add exponential sampling
- int a = 0;
- next_step = -int_length * log(s_rndm_(a));
} else
- next_step = std::numeric_limits<double>::infinity();
+ int_length = std::numeric_limits<double>::infinity();
/*
what are the units of the output? slant depth or 3space length?
*/
- std::cout << "ProcessSplit: "
- << "next interaction (g/cm2): " << next_step << std::endl;
- return next_step;
+ return int_length;
+ //
+ // int a = 0;
+ // const double next_step = -int_length * log(s_rndm_(a));
+ // std::cout << "ProcessSplit: "
+ // << "next step (g/cm2): " << next_step << std::endl;
+ // return next_step;
}
- template <typename Particle, typename Stack>
- EProcessReturn DoContinuous(Particle&, setup::Trajectory&, Stack&) const {
- // corsika::utls::ignore(p);
+ template <typename Particle, typename Track, typename Stack>
+ EProcessReturn DoContinuous(Particle&, Track&, Stack&) const {
return EProcessReturn::eOk;
}
template <typename Particle, typename Stack>
- void DoDiscrete(Particle& p, Stack& s) const {
- cout << "ProcessSplit: "
- << "DoDiscrete: " << p.GetPID() << " interaction? "
+ void DoInteraction(Particle& p, Stack& s) const {
+ cout << "ProcessSibyll: "
+ << "DoInteraction: " << p.GetPID() << " interaction? "
<< process::sibyll::CanInteract(p.GetPID()) << endl;
if (process::sibyll::CanInteract(p.GetPID())) {
cout << "defining coordinates" << endl;
@@ -351,10 +360,10 @@ public:
// get energy of particle from stack
/*
- stack is in GeV in lab. frame
- convert to GeV in cm. frame
- (assuming proton at rest as target AND
- assuming no pT, i.e. shower frame-z is aligned with hadron-int-frame-z)
+ stack is in GeV in lab. frame
+ convert to GeV in cm. frame
+ (assuming proton at rest as target AND
+ assuming no pT, i.e. shower frame-z is aligned with hadron-int-frame-z)
*/
// total energy: E_beam + E_target
// in lab. frame: E_beam + m_target*c**2
@@ -382,11 +391,13 @@ public:
int kBeam = process::sibyll::ConvertToSibyllRaw(p.GetPID());
std::cout << "ProcessSplit: "
- << " DoDiscrete: E(GeV):" << E / 1_GeV << " Ecm(GeV): " << Ecm / 1_GeV
+ << " DoInteraction: E(GeV):" << E / 1_GeV << " Ecm(GeV): " << Ecm / 1_GeV
<< std::endl;
if (E < 8.5_GeV || Ecm < 10_GeV) {
std::cout << "ProcessSplit: "
- << " DoDiscrete: low en. particle, skipping.." << std::endl;
+ << " DoInteraction: dropping particle.." << std::endl;
+ p.Delete();
+ fCount++;
} else {
// Sibyll does not know about units..
double sqs = Ecm / 1_GeV;
@@ -458,9 +469,8 @@ public:
fCount = 0;
corsika::random::RNGManager& rmng = corsika::random::RNGManager::GetInstance();
- ;
- const std::string str_name = "s_rndm";
- rmng.RegisterRandomStream(str_name);
+
+ rmng.RegisterRandomStream("s_rndm");
// test random number generator
std::cout << "ProcessSplit: "
@@ -474,7 +484,6 @@ public:
setTrackedParticlesStable();
}
-
int GetCount() { return fCount; }
EnergyType GetEnergy() { return fEnergy; }
@@ -488,12 +497,27 @@ double s_rndm_(int&) {
return rmng() / (double)rmng.max();
}
-
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));
CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
- tracking_line::TrackingLine<setup::Stack> tracking;
+ tracking_line::TrackingLine<setup::Stack> tracking(env);
stack_inspector::StackInspector<setup::Stack> p0(true);
ProcessSplit p1;
diff --git a/Environment/CMakeLists.txt b/Environment/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..420fc20f087d0cb8ba3ab884f4ae0db7efc8c469
--- /dev/null
+++ b/Environment/CMakeLists.txt
@@ -0,0 +1,52 @@
+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)
+
diff --git a/Environment/Environment.h b/Environment/Environment.h
new file mode 100644
index 0000000000000000000000000000000000000000..5afe69c57000976ef529559696e6bcdca31da609
--- /dev/null
+++ b/Environment/Environment.h
@@ -0,0 +1,57 @@
+/**
+ * (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
diff --git a/Environment/HomogeneousMedium.h b/Environment/HomogeneousMedium.h
new file mode 100644
index 0000000000000000000000000000000000000000..a8fa5f7283cfe0e776383f24915b98748d54349f
--- /dev/null
+++ b/Environment/HomogeneousMedium.h
@@ -0,0 +1,59 @@
+/**
+ * (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
diff --git a/Environment/IMediumModel.h b/Environment/IMediumModel.h
new file mode 100644
index 0000000000000000000000000000000000000000..4a1eedd86562187c310d180a4516ee3425a05674
--- /dev/null
+++ b/Environment/IMediumModel.h
@@ -0,0 +1,34 @@
+#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
diff --git a/Environment/NuclearComposition.h b/Environment/NuclearComposition.h
new file mode 100644
index 0000000000000000000000000000000000000000..46a2bf5e859f10a7f5f13a4e4312d38cb1ae6086
--- /dev/null
+++ b/Environment/NuclearComposition.h
@@ -0,0 +1,36 @@
+#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
diff --git a/Environment/VolumeTreeNode.h b/Environment/VolumeTreeNode.h
new file mode 100644
index 0000000000000000000000000000000000000000..d92d77b398bd180d9d64dc5039f32e50ca9be487
--- /dev/null
+++ b/Environment/VolumeTreeNode.h
@@ -0,0 +1,119 @@
+/**
+ * (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
diff --git a/Environment/testEnvironment.cc b/Environment/testEnvironment.cc
new file mode 100644
index 0000000000000000000000000000000000000000..8a8e9c13b5a899ad983d0e852379e5621c9f1489
--- /dev/null
+++ b/Environment/testEnvironment.cc
@@ -0,0 +1,30 @@
+/**
+ * (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);
+}
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index 4f4aadf786f7f7b8848cac9e6c2370e0ad19c5b9..7cd9d077de0bc6f84174789771714978d8ac3f94 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -13,12 +13,13 @@
#define _include_corsika_cascade_Cascade_h_
#include <corsika/process/ProcessReturn.h>
+#include <corsika/random/RNGManager.h>
+#include <corsika/setup/SetupTrajectory.h>
#include <corsika/units/PhysicalUnits.h>
#include <type_traits>
-#include <corsika/setup/SetupTrajectory.h>
-
+using namespace corsika;
using namespace corsika::units::si;
namespace corsika::cascade {
@@ -50,26 +51,92 @@ namespace corsika::cascade {
}
// do cascade equations, which can put new particles on Stack,
// thus, the double loop
- // DoCascadeEquations(); //
+ // DoCascadeEquations();
}
}
void Step(Particle& particle) {
- std::cout << "+++++++++++++++++++++++++++++++" << std::endl;
- std::cout << "Cascade: starting step.." << std::endl;
+
+ // get access to random number generator
+ static corsika::random::RNG& rmng =
+ corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
+
+ // determine geometric tracking
corsika::setup::Trajectory step = fTracking.GetTrack(particle);
- fProcesseList.MinStepLength(particle, step);
+
+ // determine combined total interaction length (inverse)
+ const double total_inv_lambda =
+ fProcesseList.GetTotalInverseInteractionLength(particle, step);
+
+ // sample random exponential step length
+ const double sample_step = rmng() / (double)rmng.max();
+ const double next_interact = -log(sample_step) / total_inv_lambda;
+ std::cout << "total_inv_lambda=" << total_inv_lambda
+ << ", next_interact=" << next_interact << std::endl;
+
+ // determine the maximum geometric step length
+ const double distance_max = fProcesseList.MaxStepLength(particle, step);
+ std::cout << "distance_max=" << distance_max << std::endl;
+
+ // determine combined total inverse decay time
+ const double total_inv_lifetime = fProcesseList.GetTotalInverseLifetime(particle);
+
+ // sample random exponential decay time
+ const double sample_decay = rmng() / (double)rmng.max();
+ const double next_decay = -log(sample_decay) / total_inv_lifetime;
+ std::cout << "total_inv_lifetime=" << total_inv_lifetime
+ << ", next_decay=" << next_decay << std::endl;
+
+ // convert next_step from grammage to length [m]
+ // Environment::GetDistance(step, next_step);
+ const double distance_interact = next_interact;
+ // ....
+
+ // convert next_decay from time to length [m]
+ // Environment::GetDistance(step, next_decay);
+ const double distance_decay = next_decay;
+ // ....
+
+ // take minimum of geometry, interaction, decay for next step
+ const double distance_decay_interact = std::min(next_decay, next_interact);
+ const double distance_next = std::min(distance_decay_interact, distance_max);
/// here the particle is actually moved along the trajectory to new position:
// std::visit(corsika::setup::ParticleUpdate<Particle>{particle}, step);
particle.SetPosition(step.GetPosition(1));
+ // .... also update time, momentum, direction, ...
+ // apply all continuous processes on particle + track
corsika::process::EProcessReturn status =
fProcesseList.DoContinuous(particle, step, fStack);
+
if (status == corsika::process::EProcessReturn::eParticleAbsorbed) {
- fStack.Delete(particle); // TODO: check if this is really needed
+ // fStack.Delete(particle); // TODO: check if this is really needed
} else {
- fProcesseList.DoDiscrete(particle, fStack);
+
+ std::cout << "select process " << (distance_decay_interact < distance_max)
+ << std::endl;
+ // check if geometry limits step, then quit this step
+ if (distance_decay_interact < distance_max) {
+ // check weather decay or interaction limits this step
+ if (distance_decay > distance_interact) {
+ std::cout << "collide" << std::endl;
+ const double actual_inv_length =
+ fProcesseList.GetTotalInverseInteractionLength(particle, step);
+ const double sample_process = rmng() / (double)rmng.max();
+ double inv_lambda_count = 0;
+ fProcesseList.SelectInteraction(particle, fStack, actual_inv_length,
+ sample_process, inv_lambda_count);
+ } else {
+ std::cout << "decay" << std::endl;
+ const double actual_decay_time =
+ fProcesseList.GetTotalInverseLifetime(particle);
+ const double sample_process = rmng() / (double)rmng.max();
+ double inv_decay_count = 0;
+ fProcesseList.SelectDecay(particle, fStack, actual_decay_time, sample_process,
+ inv_decay_count);
+ }
+ }
}
}
diff --git a/Framework/Cascade/SibStack.h b/Framework/Cascade/SibStack.h
index aed25d58c615a3f213781d92dae562783e4eecd4..a25316553da2182f6da9b2668f7f67336ede9a84 100644
--- a/Framework/Cascade/SibStack.h
+++ b/Framework/Cascade/SibStack.h
@@ -6,7 +6,6 @@
#include <corsika/cascade/sibyll2.3c.h>
#include <corsika/process/sibyll/ParticleConversion.h>
-#include <corsika/units/PhysicalUnits.h>
#include <corsika/stack/Stack.h>
#include <corsika/units/PhysicalUnits.h>
diff --git a/Framework/Cascade/testCascade.cc b/Framework/Cascade/testCascade.cc
index a218cb788e9da5748fbc634367cce8810a5356c8..423940c156782b2978ed402524f209cf489178fb 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,10 @@
* the license.
*/
+#include <limits>
+
+#include <corsika/environment/Environment.h>
+
#include <corsika/cascade/Cascade.h>
#include <corsika/process/ProcessSequence.h>
@@ -17,6 +20,8 @@
#include <corsika/stack/super_stupid/SuperStupidStack.h>
+#include <corsika/particles/ParticleProperties.h>
+
#include <corsika/geometry/Point.h>
#include <corsika/geometry/RootCoordinateSystem.h>
#include <corsika/geometry/Vector.h>
@@ -36,23 +41,21 @@ using namespace corsika::geometry;
#include <iostream>
using namespace std;
+using namespace corsika::units::si;
static int fCount = 0;
-class ProcessSplit : public corsika::process::BaseProcess<ProcessSplit> {
+class ProcessSplit : public corsika::process::ContinuousProcess<ProcessSplit> {
public:
ProcessSplit() {}
- template <typename Particle>
- void MinStepLength(Particle&, setup::Trajectory&) const {}
-
- template <typename Particle, typename Stack>
- EProcessReturn DoContinuous(Particle&, setup::Trajectory&, Stack&) const {
- return EProcessReturn::eOk;
+ template <typename Particle, typename T>
+ double MaxStepLength(Particle&, T&) const {
+ return 1;
}
- template <typename Particle, typename Stack>
- void DoDiscrete(Particle& p, Stack& s) const {
+ template <typename Particle, typename T, typename Stack>
+ void DoContinuous(Particle& p, T&, Stack& s) const {
EnergyType E = p.GetEnergy();
if (E < 85_MeV) {
p.Delete();
@@ -60,7 +63,9 @@ public:
} else {
p.SetEnergy(E / 2);
auto pnew = s.NewParticle();
- // s.Copy(p, pnew);
+ // s.Copy(p, pnew); fix that .... todo
+ pnew.SetPID(p.GetPID());
+ pnew.SetTime(p.GetTime());
pnew.SetEnergy(E / 2);
pnew.SetPosition(p.GetPosition());
pnew.SetMomentum(p.GetMomentum());
@@ -76,7 +81,20 @@ private:
TEST_CASE("Cascade", "[Cascade]") {
- tracking_line::TrackingLine<setup::Stack> tracking;
+ corsika::random::RNGManager& rmng = corsika::random::RNGManager::GetInstance();
+ ;
+ const std::string str_name = "s_rndm";
+ rmng.RegisterRandomStream(str_name);
+
+ 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;
@@ -84,19 +102,21 @@ TEST_CASE("Cascade", "[Cascade]") {
setup::Stack stack;
corsika::cascade::Cascade EAS(tracking, sequence, stack);
-
CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
stack.Clear();
auto particle = stack.NewParticle();
EnergyType E0 = 100_GeV;
+ particle.SetPID(particles::Code::Electron);
particle.SetEnergy(E0);
particle.SetPosition(Point(rootCS, {0_m, 0_m, 10_km}));
particle.SetMomentum(corsika::stack::super_stupid::MomentumVector(
rootCS, {0 * newton * second, 0 * newton * second, -1 * newton * second}));
+ particle.SetTime(0_ns);
EAS.Init();
EAS.Run();
+ /*
SECTION("sectionTwo") {
for (int i = 0; i < 0; ++i) {
stack.Clear();
@@ -109,4 +129,5 @@ TEST_CASE("Cascade", "[Cascade]") {
// cout << "Result: E0=" << E0 / 1_GeV << "GeV, count=" << p1.GetCount() << endl;
}
}
+ */
}
diff --git a/Framework/Geometry/BaseTrajectory.h b/Framework/Geometry/BaseTrajectory.h
index 8289e0dca67fbdaa2865fd7c7799ccf51d619c35..5540fc89d62e240217b79346d86aa32b04004baf 100644
--- a/Framework/Geometry/BaseTrajectory.h
+++ b/Framework/Geometry/BaseTrajectory.h
@@ -42,8 +42,11 @@ namespace corsika::geometry {
* parameterized by \arg t1 and \arg t2. Requires \arg t2 > \arg t1.
*/
- virtual LengthType GetDistance(corsika::units::si::TimeType t1,
- corsika::units::si::TimeType t2) const = 0;
+ 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(
corsika::units::si::TimeType t1, corsika::units::si::TimeType t2) const {
diff --git a/Framework/Geometry/CMakeLists.txt b/Framework/Geometry/CMakeLists.txt
index 2b0c8d3af13a5ba00b5fbe6674febb0b21461f96..3e2a8e1d4298af618835e80a9f1f1751704bbe18 100644
--- a/Framework/Geometry/CMakeLists.txt
+++ b/Framework/Geometry/CMakeLists.txt
@@ -10,6 +10,7 @@ set (
Point.h
Line.h
Sphere.h
+ Volume.h
CoordinateSystem.h
RootCoordinateSystem.h
Helix.h
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/Framework/Geometry/Helix.h b/Framework/Geometry/Helix.h
index 8eb215820097fae23eef0e535631384d31db01c4..2c8f6b8792cb204e652d1222b38f7849b70e0792 100644
--- a/Framework/Geometry/Helix.h
+++ b/Framework/Geometry/Helix.h
@@ -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
diff --git a/Framework/Geometry/Line.h b/Framework/Geometry/Line.h
index e093a0055fb82ceaeab2c656227b2c005f15958e..a79ba9ce85167ebdf5cef65622e060e844cc22b4 100644
--- a/Framework/Geometry/Line.h
+++ b/Framework/Geometry/Line.h
@@ -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
diff --git a/Framework/Geometry/Sphere.h b/Framework/Geometry/Sphere.h
index 6ab7c8b9690a7d92e51a6cbb75fc3a73fc227f8c..3e5458b7e2328d062e61307f4c01589c8ee02f05 100644
--- a/Framework/Geometry/Sphere.h
+++ b/Framework/Geometry/Sphere.h
@@ -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
diff --git a/Framework/Geometry/Trajectory.h b/Framework/Geometry/Trajectory.h
index 686799756504138e349aacfa428b325166ef1e23..5bfcea0216ee6fd496275003b307cccc3f804ee8 100644
--- a/Framework/Geometry/Trajectory.h
+++ b/Framework/Geometry/Trajectory.h
@@ -12,6 +12,7 @@
#ifndef _include_TRAJECTORY_H
#define _include_TRAJECTORY_H
+#include <corsika/geometry/Point.h>
#include <corsika/units/PhysicalUnits.h>
using corsika::units::si::LengthType;
@@ -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);
}
};
diff --git a/Framework/Geometry/Vector.h b/Framework/Geometry/Vector.h
index cb4e4900544705c10a37ff16a0fd4ba00b0f9b29..94394b90dd859761d4f15d0a2389f173726ee34d 100644
--- a/Framework/Geometry/Vector.h
+++ b/Framework/Geometry/Vector.h
@@ -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
diff --git a/Framework/Geometry/Volume.h b/Framework/Geometry/Volume.h
new file mode 100644
index 0000000000000000000000000000000000000000..ac892e50b926758d76298675a3c0db80557ab9a4
--- /dev/null
+++ b/Framework/Geometry/Volume.h
@@ -0,0 +1,19 @@
+#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
diff --git a/Framework/Geometry/testGeometry.cc b/Framework/Geometry/testGeometry.cc
index 9e3be3310bbaf3ed6c9ae22bb6faa67701f20687..8130114aacb6c0bbac07b0f42f3d043d704b6c89 100644
--- a/Framework/Geometry/testGeometry.cc
+++ b/Framework/Geometry/testGeometry.cc
@@ -132,7 +132,15 @@ TEST_CASE("Sphere") {
Point center(rootCS, {0_m, 3_m, 4_m});
Sphere sphere(center, 5_m);
- SECTION("isInside") {
+ SECTION("GetCenter") {
+ CHECK((sphere.GetCenter().GetCoordinates(rootCS) -
+ QuantityVector<length_d>(0_m, 3_m, 4_m))
+ .norm()
+ .magnitude() == Approx(0).margin(absMargin));
+ CHECK(sphere.GetRadius() / 5_m == Approx(1));
+ }
+
+ SECTION("Contains") {
REQUIRE_FALSE(sphere.Contains(Point(rootCS, {100_m, 0_m, 0_m})));
REQUIRE(sphere.Contains(Point(rootCS, {2_m, 3_m, 4_m})));
}
@@ -152,11 +160,11 @@ TEST_CASE("Trajectories") {
.norm()
.magnitude() == Approx(0).margin(absMargin));
- Trajectory<Line> base(line, 1_s);
- CHECK(line.GetPosition(2_s).GetCoordinates() ==
- base.GetPosition(2_s).GetCoordinates());
+ auto const t = 1_s;
+ Trajectory<Line> base(line, t);
+ CHECK(line.GetPosition(t).GetCoordinates() == base.GetPosition(1.).GetCoordinates());
- CHECK(base.GetDistanceBetween(1_s, 2_s) / 1_m == Approx(1));
+ CHECK(base.ArcLength(1_s, 2_s) / 1_m == Approx(1));
}
SECTION("Helix") {
@@ -180,10 +188,10 @@ TEST_CASE("Trajectories") {
.norm()
.magnitude() == Approx(0).margin(absMargin));
- Trajectory<Helix> const base(helix, 1_s);
- CHECK(helix.GetPosition(1234_s).GetCoordinates() ==
- base.GetPosition(1234_s).GetCoordinates());
+ auto const t = 1234_s;
+ Trajectory<Helix> const base(helix, t);
+ CHECK(helix.GetPosition(t).GetCoordinates() == base.GetPosition(1.).GetCoordinates());
- CHECK(base.GetDistanceBetween(1_s, 2_s) / 1_m == Approx(5));
+ CHECK(base.ArcLength(1_s, 2_s) / 1_m == Approx(5));
}
}
diff --git a/Framework/ProcessSequence/BaseProcess.h b/Framework/ProcessSequence/BaseProcess.h
index 34dcb198e4d26f409837d253f3270cc0f488098b..9cc29196be0d136598e2f32480acd35dd16b158e 100644
--- a/Framework/ProcessSequence/BaseProcess.h
+++ b/Framework/ProcessSequence/BaseProcess.h
@@ -29,14 +29,29 @@ namespace corsika::process {
struct BaseProcess {
derived& GetRef() { return static_cast<derived&>(*this); }
const derived& GetRef() const { return static_cast<const derived&>(*this); }
-
+ /*
template <typename Particle, typename Stack>
inline EProcessReturn DoDiscrete(Particle&, Stack&) const; // {}
template <typename Particle, typename Track, typename Stack>
inline EProcessReturn DoContinuous(Particle&, Track&, Stack&) const; // {}
+ */
+ /*
+ template <typename Particle, typename Track>
+ inline double GetInverseInteractionLength(Particle& p, Track& t) const {
+ return 1./GetRef().GetInteractionLength(p, t);
+ }
+ */
};
+ /*
+ template<template<typename, typename> class T, typename A, typename B>
+ typename BaseProcess< T<A, B> >::is_process_sequence
+ {
+ static const bool value = true;
+ };
+ */
+
/*
template <typename T>
struct is_base {
diff --git a/Framework/ProcessSequence/CMakeLists.txt b/Framework/ProcessSequence/CMakeLists.txt
index c1f5a45862faffc92aba606e8fc1abde94220d89..6c3e5ff8cc26e571955b937920dc086fd925869c 100644
--- a/Framework/ProcessSequence/CMakeLists.txt
+++ b/Framework/ProcessSequence/CMakeLists.txt
@@ -12,7 +12,8 @@ set (
CORSIKAprocesssequence_HEADERS
BaseProcess.h
ContinuousProcess.h
- DiscreteProcess.h
+ InteractionProcess.h
+ DecayProcess.h
ProcessSequence.h
ProcessReturn.h
)
diff --git a/Framework/ProcessSequence/DecayProcess.h b/Framework/ProcessSequence/DecayProcess.h
new file mode 100644
index 0000000000000000000000000000000000000000..8a49803ed49913a1c8cf5c217a5c866de307f4bd
--- /dev/null
+++ b/Framework/ProcessSequence/DecayProcess.h
@@ -0,0 +1,51 @@
+
+/**
+ * (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_corsika_decayprocess_h_
+#define _include_corsika_decayprocess_h_
+
+#include <corsika/process/ProcessReturn.h> // for convenience
+#include <corsika/setup/SetupTrajectory.h>
+
+namespace corsika::process {
+
+ /**
+ \class DecayProcess
+
+ The structural base type of a process object in a
+ ProcessSequence. Both, the ProcessSequence and all its elements
+ are of type DecayProcess<T>
+
+ */
+
+ template <typename derived>
+ struct DecayProcess {
+
+ derived& GetRef() { return static_cast<derived&>(*this); }
+ const derived& GetRef() const { return static_cast<const derived&>(*this); }
+
+ /// here starts the interface-definition part
+ // -> enforce derived to implement DoDecay...
+ template <typename Particle, typename Stack>
+ inline EProcessReturn DoDecay(Particle&, Stack&) const;
+
+ template <typename Particle>
+ inline double GetLifetime(Particle& p) const;
+
+ template <typename Particle>
+ inline double GetInverseLifetime(Particle& p) const {
+ return 1. / GetRef().GetLifetime(p);
+ }
+ };
+
+} // namespace corsika::process
+
+#endif
diff --git a/Framework/ProcessSequence/DiscreteProcess.h b/Framework/ProcessSequence/DiscreteProcess.h
index a540c92c783dcad931f0d2bf8278c057736bda79..7ce50810f7fccbbbcb0c707d4124ec1ad03f0244 100644
--- a/Framework/ProcessSequence/DiscreteProcess.h
+++ b/Framework/ProcessSequence/DiscreteProcess.h
@@ -36,7 +36,15 @@ namespace corsika::process {
/// here starts the interface-definition part
// -> enforce derived to implement DoDiscrete...
template <typename Particle, typename Stack>
- inline EProcessReturn DoDiscrete(Particle&, Stack&) const; // {}
+ inline EProcessReturn DoDiscrete(Particle&, Stack&) const;
+
+ template <typename Particle, typename Track>
+ inline double GetInteractionLength(Particle& p, Track& t) const;
+
+ template <typename Particle, typename Track>
+ inline double GetInverseInteractionLength(Particle& p, Track& t) const {
+ return 1. / GetRef().GetInteractionLength(p, t);
+ }
};
} // namespace corsika::process
diff --git a/Framework/ProcessSequence/InteractionProcess.h b/Framework/ProcessSequence/InteractionProcess.h
new file mode 100644
index 0000000000000000000000000000000000000000..40264e606a7de8ecf607fa8d5173190e9e4ca0e4
--- /dev/null
+++ b/Framework/ProcessSequence/InteractionProcess.h
@@ -0,0 +1,51 @@
+
+/**
+ * (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_corsika_interactionprocess_h_
+#define _include_corsika_interactionprocess_h_
+
+#include <corsika/process/ProcessReturn.h> // for convenience
+#include <corsika/setup/SetupTrajectory.h>
+
+namespace corsika::process {
+
+ /**
+ \class InteractionProcess
+
+ The structural base type of a process object in a
+ ProcessSequence. Both, the ProcessSequence and all its elements
+ are of type InteractionProcess<T>
+
+ */
+
+ template <typename derived>
+ struct InteractionProcess {
+
+ derived& GetRef() { return static_cast<derived&>(*this); }
+ const derived& GetRef() const { return static_cast<const derived&>(*this); }
+
+ /// here starts the interface-definition part
+ // -> enforce derived to implement DoInteraction...
+ template <typename Particle, typename Stack>
+ inline EProcessReturn DoInteraction(Particle&, Stack&) const;
+
+ template <typename Particle, typename Track>
+ inline double GetInteractionLength(Particle& p, Track& t) const;
+
+ template <typename Particle, typename Track>
+ inline double GetInverseInteractionLength(Particle& p, Track& t) const {
+ return 1. / GetRef().GetInteractionLength(p, t);
+ }
+ };
+
+} // namespace corsika::process
+
+#endif
diff --git a/Framework/ProcessSequence/ProcessReturn.h b/Framework/ProcessSequence/ProcessReturn.h
index 82cc816c5ae9bd9d731f00c33f9bf0321a02d650..afd75d690a3b5c1ee54d7e0884a536f42221d2a2 100644
--- a/Framework/ProcessSequence/ProcessReturn.h
+++ b/Framework/ProcessSequence/ProcessReturn.h
@@ -23,6 +23,8 @@ namespace corsika::process {
enum class EProcessReturn {
eOk = 1,
eParticleAbsorbed = 2,
+ eInteracted = 3,
+ eDecayed = 4,
};
} // namespace corsika::process
diff --git a/Framework/ProcessSequence/ProcessSequence.h b/Framework/ProcessSequence/ProcessSequence.h
index 54ac90e7a18f7c40c1e5a53c82ebf854e5abe486..8517d4345202f8fb667f830eab3c093937b2ad05 100644
--- a/Framework/ProcessSequence/ProcessSequence.h
+++ b/Framework/ProcessSequence/ProcessSequence.h
@@ -14,9 +14,14 @@
#include <corsika/process/BaseProcess.h>
#include <corsika/process/ContinuousProcess.h>
-#include <corsika/process/DiscreteProcess.h>
+#include <corsika/process/DecayProcess.h>
+//#include <corsika/process/DiscreteProcess.h>
+#include <corsika/process/InteractionProcess.h>
#include <corsika/process/ProcessReturn.h>
+#include <cmath>
+#include <limits>
+
//#include <corsika/setup/SetupTrajectory.h>
// using corsika::setup::Trajectory;
//#include <variant>
@@ -25,7 +30,7 @@
namespace corsika::process {
- /* namespace detail { */
+ // namespace detail {
/* /\* template<typename TT1, typename TT2, typename Type = void> *\/ */
/* /\* struct CallHello { *\/ */
@@ -41,7 +46,7 @@ namespace corsika::process {
/* /\* static void Call(const TT1&, const TT2&) { *\/ */
/* /\* std::cout << "special" << std::endl; *\/ */
/* /\* } *\/ */
- /* /\* }; *\/ */
+ /* }; */
/* template<typename T1, typename T2, typename Particle, typename Trajectory, typename
* Stack> //, typename Type = void> */
@@ -121,7 +126,7 @@ namespace corsika::process {
/* } */
/* }; */
/* *\/ */
- /* } // end namespace detail */
+ //} // end namespace detail
/**
\class ProcessSequence
@@ -134,8 +139,15 @@ namespace corsika::process {
https://en.wikipedia.org/wiki/Curiously_recurring_template_pattern
*/
+ // this is a marker to track which BaseProcess is also a ProcessSequence
+ template <typename T>
+ struct is_process_sequence {
+ static const bool value = false;
+ };
+
template <typename T1, typename T2>
class ProcessSequence : public BaseProcess<ProcessSequence<T1, T2> > {
+
public:
const T1& A;
const T2& B;
@@ -150,41 +162,154 @@ namespace corsika::process {
template <typename Particle, typename Track, typename Stack>
inline EProcessReturn DoContinuous(Particle& p, Track& t, Stack& s) const {
EProcessReturn ret = EProcessReturn::eOk;
- if constexpr (!std::is_base_of<DiscreteProcess<T1>, T1>::value) {
- // A.DoContinuous(std::forward<Particle>(p), t, std::forward<Stack>(s));
+ if constexpr (std::is_base_of<ContinuousProcess<T1>, T1>::value ||
+ is_process_sequence<T1>::value) {
A.DoContinuous(p, t, s);
}
- if constexpr (!std::is_base_of<DiscreteProcess<T2>, T2>::value) {
- // B.DoContinuous(std::forward<Particle>(p), t, std::forward<Stack>(s));
+ if constexpr (std::is_base_of<ContinuousProcess<T2>, T2>::value ||
+ is_process_sequence<T2>::value) {
B.DoContinuous(p, t, s);
}
return ret;
}
template <typename Particle, typename Track>
- inline void MinStepLength(Particle& p, Track& track) const {
- A.MinStepLength(p, track);
- B.MinStepLength(p, track);
+ inline double MaxStepLength(Particle& p, Track& track) const {
+ double max_length = std::numeric_limits<double>::infinity();
+ if constexpr (std::is_base_of<ContinuousProcess<T1>, T1>::value ||
+ is_process_sequence<T1>::value) {
+ max_length = std::min(max_length, A.MaxStepLength(p, track));
+ }
+ if constexpr (std::is_base_of<ContinuousProcess<T2>, T2>::value ||
+ is_process_sequence<T2>::value) {
+ max_length = std::min(max_length, B.MaxStepLength(p, track));
+ }
+ return max_length;
}
- /*
template <typename Particle, typename Track>
- inline Track Transport(Particle& p, double& length) const {
- A.Transport(p, length); // todo: maybe check (?) if there is more than one Transport
- // process implemented??
- return B.Transport(
- p, length); // need to do this also to decide which Track to return!!!!
+ inline double GetTotalInteractionLength(Particle& p, Track& t) const {
+ return 1. / GetInverseInteractionLength(p, t);
+ }
+
+ template <typename Particle, typename Track>
+ inline double GetTotalInverseInteractionLength(Particle& p, Track& t) const {
+ return GetInverseInteractionLength(p, t);
+ }
+
+ template <typename Particle, typename Track>
+ inline double GetInverseInteractionLength(Particle& p, Track& t) const {
+ double tot = 0;
+ if constexpr (std::is_base_of<InteractionProcess<T1>, T1>::value ||
+ is_process_sequence<T1>::value) {
+ tot += A.GetInverseInteractionLength(p, t);
+ }
+ if constexpr (std::is_base_of<InteractionProcess<T2>, T2>::value ||
+ is_process_sequence<T2>::value) {
+ tot += B.GetInverseInteractionLength(p, t);
+ }
+ return tot;
}
- */
template <typename Particle, typename Stack>
- inline EProcessReturn DoDiscrete(Particle& p, Stack& s) const {
- if constexpr (!std::is_base_of<ContinuousProcess<T1>, T1>::value) {
- A.DoDiscrete(p, s);
+ inline EProcessReturn SelectInteraction(Particle& p, Stack& s,
+ const double lambda_inv_tot,
+ const double rndm_select,
+ double& lambda_inv_count) const {
+ if constexpr (is_process_sequence<T1>::value) {
+ // if A is a process sequence --> check inside
+ const EProcessReturn ret =
+ A.SelectInteraction(p, s, lambda_inv_count, rndm_select, lambda_inv_count);
+ // if A did succeed, stop routine
+ if (ret != EProcessReturn::eOk) { return ret; }
+ } else if constexpr (std::is_base_of<InteractionProcess<T1>, T1>::value) {
+ // if this is not a ContinuousProcess --> evaluate probability
+ lambda_inv_count += A.GetInverseInteractionLength(p, s);
+ // check if we should execute THIS process and then EXIT
+ if (rndm_select < lambda_inv_count / lambda_inv_tot) {
+ A.DoInteraction(p, s);
+ return EProcessReturn::eInteracted;
+ }
+ } // end branch A
+
+ if constexpr (is_process_sequence<T2>::value) {
+ // if A is a process sequence --> check inside
+ const EProcessReturn ret =
+ B.SelectInteraction(p, s, lambda_inv_count, rndm_select, lambda_inv_count);
+ // if A did succeed, stop routine
+ if (ret != EProcessReturn::eOk) { return ret; }
+ } else if constexpr (std::is_base_of<InteractionProcess<T2>, T2>::value) {
+ // if this is not a ContinuousProcess --> evaluate probability
+ lambda_inv_count += B.GetInverseInteractionLength(p, s);
+ // check if we should execute THIS process and then EXIT
+ if (rndm_select < lambda_inv_count / lambda_inv_tot) {
+ B.DoInteraction(p, s);
+ return EProcessReturn::eInteracted;
+ }
+ } // end branch A
+ return EProcessReturn::eOk;
+ }
+
+ template <typename Particle>
+ inline double GetTotalLifetime(Particle& p) const {
+ return 1. / GetInverseLifetime(p);
+ }
+
+ template <typename Particle>
+ inline double GetTotalInverseLifetime(Particle& p) const {
+ return GetInverseLifetime(p);
+ }
+
+ template <typename Particle>
+ inline double GetInverseLifetime(Particle& p) const {
+ double tot = 0;
+ if constexpr (std::is_base_of<DecayProcess<T1>, T1>::value ||
+ is_process_sequence<T1>::value) {
+ tot += A.GetInverseLifetime(p);
}
- if constexpr (!std::is_base_of<ContinuousProcess<T2>, T2>::value) {
- B.DoDiscrete(p, s);
+ if constexpr (std::is_base_of<DecayProcess<T2>, T2>::value ||
+ is_process_sequence<T2>::value) {
+ tot += B.GetInverseLifetime(p);
}
+ return tot;
+ }
+
+ // select decay process
+ template <typename Particle, typename Stack>
+ inline EProcessReturn SelectDecay(Particle& p, Stack& s, const double decay_inv_tot,
+ const double rndm_select,
+ double& decay_inv_count) const {
+ if constexpr (is_process_sequence<T1>::value) {
+ // if A is a process sequence --> check inside
+ const EProcessReturn ret =
+ A.SelectDecay(p, s, decay_inv_count, rndm_select, decay_inv_count);
+ // if A did succeed, stop routine
+ if (ret != EProcessReturn::eOk) { return ret; }
+ } else if constexpr (std::is_base_of<DecayProcess<T1>, T1>::value) {
+ // if this is not a ContinuousProcess --> evaluate probability
+ decay_inv_count += A.GetInverseLifetime(p);
+ // check if we should execute THIS process and then EXIT
+ if (rndm_select < decay_inv_count / decay_inv_tot) {
+ A.DoDecay(p, s);
+ return EProcessReturn::eDecayed;
+ }
+ } // end branch A
+
+ if constexpr (is_process_sequence<T2>::value) {
+ // if A is a process sequence --> check inside
+ const EProcessReturn ret =
+ B.SelectDecay(p, s, decay_inv_count, rndm_select, decay_inv_count);
+ // if A did succeed, stop routine
+ if (ret != EProcessReturn::eOk) { return ret; }
+ } else if constexpr (std::is_base_of<DecayProcess<T2>, T2>::value) {
+ // if this is not a ContinuousProcess --> evaluate probability
+ decay_inv_count += B.GetInverseLifetime(p);
+ // check if we should execute THIS process and then EXIT
+ if (rndm_select < decay_inv_count / decay_inv_tot) {
+ B.DoDecay(p, s);
+ return EProcessReturn::eDecayed;
+ }
+ } // end branch B
return EProcessReturn::eOk;
}
@@ -195,8 +320,8 @@ namespace corsika::process {
}
};
- /// the +operator assembles many BaseProcess, ContinuousProcess, and
- /// DiscreteProcess objects into a ProcessSequence, all combinatorics
+ /// the + operator assembles many BaseProcess, ContinuousProcess, and
+ /// InteractionProcess objects into a ProcessSequence, all combinatorics
/// must be allowed, this is why we define a macro to define all
/// combinations here:
@@ -207,87 +332,26 @@ namespace corsika::process {
}
OPSEQ(BaseProcess, BaseProcess)
- OPSEQ(BaseProcess, DiscreteProcess)
+ OPSEQ(BaseProcess, InteractionProcess)
OPSEQ(BaseProcess, ContinuousProcess)
+ OPSEQ(BaseProcess, DecayProcess)
OPSEQ(ContinuousProcess, BaseProcess)
- OPSEQ(ContinuousProcess, DiscreteProcess)
+ OPSEQ(ContinuousProcess, InteractionProcess)
OPSEQ(ContinuousProcess, ContinuousProcess)
- OPSEQ(DiscreteProcess, BaseProcess)
- OPSEQ(DiscreteProcess, DiscreteProcess)
- OPSEQ(DiscreteProcess, ContinuousProcess)
-
- /*
- template <typename T1>
- struct depth_lhs
- {
- static const int num = 0;
- };
-
-
-
- // terminating condition
- template <typename T1, typename T2>
- struct depth_lhs< Sequence<T1,T2> >
- {
- // try to expand the left node (T1) which might be a Sequence type
- static const int num = 1 + depth_lhs<T1>::num;
- };
- */
-
- /*
- template <typename T1>
- struct mat_ptrs
- {
- static const int num = 0;
-
- inline static void
- get_ptrs(const Process** ptrs, const T1& X)
- {
- ptrs[0] = reinterpret_cast<const Process*>(&X);
- }
- };
-
-
- template <typename T1, typename T2>
- struct mat_ptrs< Sequence<T1,T2> >
- {
- static const int num = 1 + mat_ptrs<T1>::num;
-
- inline static void
- get_ptrs(const Process** in_ptrs, const Sequence<T1,T2>& X)
- {
- // traverse the left node
- mat_ptrs<T1>::get_ptrs(in_ptrs, X.A);
- // get address of the matrix on the right node
- in_ptrs[num] = reinterpret_cast<const Process*>(&X.B);
- }
- };
- */
-
- /*
- template<typename T1, typename T2>
- const Process&
- Process::operator=(const Sequence<T1,T2>& X)
- {
- int N = 1 + depth_lhs< Sequence<T1,T2> >::num;
- const Process* ptrs[N];
- mat_ptrs< Sequence<T1,T2> >::get_ptrs(ptrs, X);
- int r = ptrs[0]->rows;
- int c = ptrs[0]->cols;
- // ... check that all matrices have the same size ...
- set_size(r, c);
- for(int j=0; j<r*c; ++j)
- {
- double sum = ptrs[0]->data[j];
- for(int i=1; i<N; ++i)
- {
- sum += ptrs[i]->data[j];
- }
- data[j] = sum;
- }
- return *this;
- }
- */
+ OPSEQ(ContinuousProcess, DecayProcess)
+ OPSEQ(InteractionProcess, BaseProcess)
+ OPSEQ(InteractionProcess, InteractionProcess)
+ OPSEQ(InteractionProcess, ContinuousProcess)
+ OPSEQ(InteractionProcess, DecayProcess)
+ OPSEQ(DecayProcess, BaseProcess)
+ OPSEQ(DecayProcess, InteractionProcess)
+ OPSEQ(DecayProcess, ContinuousProcess)
+ OPSEQ(DecayProcess, DecayProcess)
+
+ template <template <typename, typename> class T, typename A, typename B>
+ struct is_process_sequence<T<A, B> > {
+ static const bool value = true;
+ };
} // namespace corsika::process
diff --git a/Framework/ProcessSequence/testProcessSequence.cc b/Framework/ProcessSequence/testProcessSequence.cc
index 7094cacb94e8da03bcfce1b9190160fb1b8c1898..af9bfeeade6c2d149f88b803f0ddd2ebb105e6aa 100644
--- a/Framework/ProcessSequence/testProcessSequence.cc
+++ b/Framework/ProcessSequence/testProcessSequence.cc
@@ -66,7 +66,7 @@ public:
}
};
-class Process1 : public DiscreteProcess<Process1> {
+class Process1 : public InteractionProcess<Process1> {
public:
Process1(const int v)
: fV(v) {}
@@ -76,7 +76,7 @@ public:
globalCount++;
}
template <typename D, typename S>
- inline EProcessReturn DoDiscrete(D& d, S&) const {
+ inline EProcessReturn DoInteraction(D& d, S&) const {
for (int i = 0; i < nData; ++i) d.p[i] += 1 + i;
return EProcessReturn::eOk;
}
@@ -84,7 +84,7 @@ public:
int fV;
};
-class Process2 : public DiscreteProcess<Process2> {
+class Process2 : public InteractionProcess<Process2> {
int fV = 0;
public:
@@ -96,13 +96,18 @@ public:
globalCount++;
}
template <typename Particle, typename Stack>
- inline EProcessReturn DoDiscrete(Particle&, Stack&) const {
- cout << "Process2::DoDiscrete" << endl;
+ inline EProcessReturn DoInteraction(Particle&, Stack&) const {
+ cout << "Process2::DoInteraction" << endl;
return EProcessReturn::eOk;
}
+ template <typename Particle, typename Track>
+ inline double GetInteractionLength(Particle&, Track&) const {
+ cout << "Process2::GetInteractionLength" << endl;
+ return 3;
+ }
};
-class Process3 : public DiscreteProcess<Process3> {
+class Process3 : public InteractionProcess<Process3> {
int fV = 0;
public:
@@ -114,10 +119,15 @@ public:
globalCount++;
}
template <typename Particle, typename Stack>
- inline EProcessReturn DoDiscrete(Particle&, Stack&) const {
- cout << "Process3::DoDiscrete" << endl;
+ inline EProcessReturn DoInteraction(Particle&, Stack&) const {
+ cout << "Process3::DoInteraction" << endl;
return EProcessReturn::eOk;
}
+ template <typename Particle, typename Track>
+ inline double GetInteractionLength(Particle&, Track&) const {
+ cout << "Process3::GetInteractionLength" << endl;
+ return 1.;
+ }
};
class Process4 : public BaseProcess<Process4> {
@@ -138,7 +148,28 @@ public:
}
// inline double MinStepLength(D& d) {
template <typename Particle, typename Stack>
- EProcessReturn DoDiscrete(Particle&, Stack&) const {
+ EProcessReturn DoInteraction(Particle&, Stack&) const {
+ return EProcessReturn::eOk;
+ }
+};
+
+class Decay1 : public DecayProcess<Decay1> {
+ int fV = 0;
+
+public:
+ Decay1(const int v)
+ : fV(v) {}
+ void Init() {
+ cout << "Decay1::Init" << endl;
+ assert(globalCount == fV);
+ globalCount++;
+ }
+ template <typename Particle>
+ double GetLifetime(Particle&) const {
+ return 1;
+ }
+ template <typename Particle, typename Stack>
+ EProcessReturn DoDecay(Particle&, Stack&) const {
return EProcessReturn::eOk;
}
};
@@ -169,6 +200,34 @@ TEST_CASE("Process Sequence", "[Process Sequence]") {
// REQUIRE_THROWS(sequence_wrong.Init());
}
+ SECTION("interaction length") {
+ ContinuousProcess1 cp1(0);
+ Process2 m2(1);
+ Process3 m3(2);
+
+ DummyStack s;
+ DummyTrajectory t;
+
+ const auto sequence2 = cp1 + m2 + m3;
+ double tot = sequence2.GetTotalInteractionLength(s, t);
+ double tot_inv = sequence2.GetTotalInverseInteractionLength(s, t);
+ cout << "lambda_tot=" << tot << " lambda_tot_inv=" << tot_inv << endl;
+ }
+
+ SECTION("lifetime") {
+ ContinuousProcess1 cp1(0);
+ Process2 m2(1);
+ Process3 m3(2);
+ Decay1 d3(2);
+
+ DummyStack s;
+
+ const auto sequence2 = cp1 + m2 + m3 + d3;
+ double tot = sequence2.GetTotalLifetime(s);
+ double tot_inv = sequence2.GetTotalInverseLifetime(s);
+ cout << "lambda_tot=" << tot << " lambda_tot_inv=" << tot_inv << endl;
+ }
+
SECTION("sectionTwo") {
ContinuousProcess1 cp1(0);
@@ -189,14 +248,14 @@ TEST_CASE("Process Sequence", "[Process Sequence]") {
sequence2.DoContinuous(p, t, s);
cout << "-->dodisc" << endl;
- sequence2.DoDiscrete(p, s);
+ // sequence2.DoInteraction(p, s);
cout << "-->done" << endl;
const int nLoop = 5;
cout << "Running loop with n=" << nLoop << endl;
for (int i = 0; i < nLoop; ++i) {
sequence2.DoContinuous(p, t, s);
- sequence2.DoDiscrete(p, s);
+ // sequence2.DoInteraction(p, s);
}
for (int i = 0; i < nData; i++) { cout << "data[" << i << "]=" << p.p[i] << endl; }
cout << "done" << endl;
diff --git a/Framework/StackInterface/Stack.h b/Framework/StackInterface/Stack.h
index b4660144df60754df6ff8e78e45d0920f5a58cf8..f47c15d3daf3071a8c71da1b81b71e575a620784 100644
--- a/Framework/StackInterface/Stack.h
+++ b/Framework/StackInterface/Stack.h
@@ -67,6 +67,7 @@ namespace corsika::stack {
IncrementSize();
return StackIterator(*this, GetSize() - 1);
}
+ void Copy(StackIterator& a, StackIterator& b) { Copy(a.GetIndex(), b.GetIndex()); }
/// delete this particle
void Delete(StackIterator& p) {
if (GetSize() == 0) { /*error*/
diff --git a/Framework/Units/PhysicalUnits.h b/Framework/Units/PhysicalUnits.h
index e6198cffd0bbdf3c443aa79b0874a82f11d318b8..7083ce743811176e2c807f2d421dde2273fbb673 100644
--- a/Framework/Units/PhysicalUnits.h
+++ b/Framework/Units/PhysicalUnits.h
@@ -52,8 +52,9 @@ namespace corsika::units::si {
using EnergyType = phys::units::quantity<phys::units::energy_d, double>;
using MassType = phys::units::quantity<phys::units::mass_d, double>;
using MassDensityType = phys::units::quantity<phys::units::mass_density_d, double>;
- using CrossSectionType = phys::units::quantity<sigma_d, double>;
+ using GrammageType = phys::units::quantity<phys::units::dimensions<-2, 1, 0>, double>;
using MomentumType = phys::units::quantity<momentum_d, double>;
+ using CrossSectionType = phys::units::quantity<sigma_d, double>;
} // end namespace corsika::units::si
diff --git a/Processes/Sibyll/ProcessDecay.h b/Processes/Sibyll/ProcessDecay.h
index e828ca9aede293be2a1787e335af661a35d01a16..701094e733ad463f673c8d73cdee0bdf34bdddab 100644
--- a/Processes/Sibyll/ProcessDecay.h
+++ b/Processes/Sibyll/ProcessDecay.h
@@ -95,7 +95,7 @@ namespace corsika::process {
const double gamma = E / m;
TimeType t0 = GetLifetime(p.GetPID());
- cout << "ProcessDecay: code: " << (p.GetPID()) << endl;
+ cout << "ProcessDecay: code: " << (p.GetPID()) << endl;
cout << "ProcessDecay: MinStep: t0: " << t0 << endl;
cout << "ProcessDecay: MinStep: gamma: " << gamma << endl;
cout << "ProcessDecay: MinStep: density: " << density << endl;
diff --git a/Processes/StackInspector/StackInspector.cc b/Processes/StackInspector/StackInspector.cc
index 949cfd025142bfb2afab60bdc7087cc9ff8a8d05..ad90442848976c4cc779977a2c2f1890cedb1091 100644
--- a/Processes/StackInspector/StackInspector.cc
+++ b/Processes/StackInspector/StackInspector.cc
@@ -10,6 +10,7 @@
*/
#include <corsika/geometry/RootCoordinateSystem.h>
+#include <corsika/particles/ParticleProperties.h>
#include <corsika/process/stack_inspector/StackInspector.h>
#include <corsika/units/PhysicalUnits.h>
@@ -18,6 +19,7 @@
#include <corsika/setup/SetupTrajectory.h>
#include <iostream>
+#include <limits>
using namespace std;
using namespace corsika;
@@ -26,7 +28,8 @@ using namespace corsika::process::stack_inspector;
template <typename Stack>
StackInspector<Stack>::StackInspector(const bool aReport)
- : fReport(aReport), fCountStep(0) {}
+ : fReport(aReport)
+ , fCountStep(0) {}
template <typename Stack>
StackInspector<Stack>::~StackInspector() {}
@@ -55,8 +58,8 @@ process::EProcessReturn StackInspector<Stack>::DoContinuous(Particle&, setup::Tr
}
template <typename Stack>
-void StackInspector<Stack>::MinStepLength(Particle&, setup::Trajectory&) const {
- // return 0;
+double StackInspector<Stack>::MaxStepLength(Particle&, setup::Trajectory&) const {
+ return std::numeric_limits<double>::infinity();
}
template <typename Stack>
diff --git a/Processes/StackInspector/StackInspector.h b/Processes/StackInspector/StackInspector.h
index dddb7e4f6fc9e1ddbfc4297b26e77e8c93a642c7..6032beac618b0941b01ffd76962c3bc13ad2300b 100644
--- a/Processes/StackInspector/StackInspector.h
+++ b/Processes/StackInspector/StackInspector.h
@@ -36,7 +36,7 @@ namespace corsika::process {
EProcessReturn DoContinuous(Particle&, corsika::setup::Trajectory&, Stack& s) const;
// template <typename Particle>
- void MinStepLength(Particle&, corsika::setup::Trajectory&) const;
+ double MaxStepLength(Particle&, corsika::setup::Trajectory&) const;
private:
bool fReport;
diff --git a/Processes/TrackingLine/CMakeLists.txt b/Processes/TrackingLine/CMakeLists.txt
index 2911c9b787787f43ba837e64fd67c89921c699b0..83ed893c150d6d041182bbe7ab9cd30046edb640 100644
--- a/Processes/TrackingLine/CMakeLists.txt
+++ b/Processes/TrackingLine/CMakeLists.txt
@@ -22,17 +22,16 @@ target_include_directories (
install (FILES ${MODEL_HEADERS} DESTINATION include/${MODEL_NAMESPACE})
-
-
-# # --------------------
-# # code unit testing
-# add_executable (testStackInspector testStackInspector.cc)
-
-# target_link_libraries (
-# testStackInspector
-# CORSIKAunits
-# CORSIKAthirdparty # for catch2
-# )
-
-# add_test (NAME testStackInspector COMMAND testStackInspector)
-
+# #-- -- -- -- -- -- -- -- -- --
+# #code unit testing
+add_executable (testTrackingLine testTrackingLine.cc)
+
+target_link_libraries (
+ testTrackingLine
+ CORSIKAunits
+ CORSIKAenvironment
+ CORSIKAgeometry
+ CORSIKAthirdparty # for catch2
+)
+
+add_test (NAME testTrackingLine COMMAND testTrackingLine)
diff --git a/Processes/TrackingLine/TrackingLine.h b/Processes/TrackingLine/TrackingLine.h
index e750342a1add1bff5e32bd270710e995d60ce4b0..b6b8864087b53fe577029c379996ca9d6ef29d35 100644
--- a/Processes/TrackingLine/TrackingLine.h
+++ b/Processes/TrackingLine/TrackingLine.h
@@ -1,14 +1,32 @@
-#ifndef _include_corsika_processes_TrackinLine_h_
-#define _include_corsika_processes_TrackinLine_h_
+/**
+ * (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_corsika_processes_TrackingLine_h_
+#define _include_corsika_processes_TrackingLine_h_
#include <corsika/geometry/Point.h>
+#include <corsika/geometry/Sphere.h>
#include <corsika/geometry/Vector.h>
#include <corsika/units/PhysicalUnits.h>
+#include <corsika/environment/Environment.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
+#include <algorithm>
+#include <iostream>
+#include <optional>
+#include <stdexcept>
+#include <utility>
+
using namespace corsika;
namespace corsika::process {
@@ -16,15 +34,96 @@ namespace corsika::process {
namespace tracking_line {
template <typename Stack>
- class TrackingLine { // Newton-step, naja.. not yet
+ class TrackingLine { //
typedef typename Stack::ParticleType Particle;
+ 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) {
+ using namespace corsika::units::si;
+ auto const& cs = fEnvironment.GetCoordinateSystem();
+ geometry::Point const origin(cs, 0_m, 0_m, 0_m);
+
+ auto const r0 = (traj.GetR0() - origin);
+ auto const v0 = traj.GetV0();
+ auto const c0 = (sphere.GetCenter() - origin);
+
+ auto const alpha = r0.dot(v0) - 2 * v0.dot(c0);
+ auto const beta = c0.squaredNorm() + r0.squaredNorm() + 2 * c0.dot(r0) -
+ sphere.GetRadius() * sphere.GetRadius();
+
+ auto const discriminant = alpha * alpha - 4 * beta * v0.squaredNorm();
+
+ //~ std::cout << "discriminant: " << discriminant << std::endl;
+ //~ std::cout << "alpha: " << alpha << std::endl;
+ //~ std::cout << "beta: " << beta << std::endl;
+
+ 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()));
+ } else {
+ return {};
+ }
+ }
+
+ TrackingLine(corsika::environment::Environment const& pEnv)
+ : fEnvironment(pEnv) {}
void Init() {}
- setup::Trajectory GetTrack(Particle& p) {
- geometry::Vector<SpeedType::dimension_type> v = p.GetDirection();
- geometry::Line traj(p.GetPosition(), v);
- return geometry::Trajectory<corsika::geometry::Line>(traj, 100_ns);
+
+ auto GetTrack(Particle& p) {
+ using namespace corsika::units::si;
+ geometry::Vector<SpeedType::dimension_type> const velocity =
+ p.GetMomentum() / p.GetEnergy() * corsika::units::si::constants::cSquared;
+
+ auto const currentPosition = p.GetPosition();
+ geometry::Line line(currentPosition, velocity);
+
+ auto const* currentVolumeNode =
+ fEnvironment.GetUniverse()->GetContainingNode(currentPosition);
+ auto const& children = currentVolumeNode->GetChildNodes();
+ auto const& excluded = currentVolumeNode->GetExcludedNodes();
+
+ std::vector<TimeType> intersectionTimes;
+
+ auto addIfIntersects = [&](auto& vtn) {
+ auto const& volume = vtn.GetVolume();
+ auto const& sphere = dynamic_cast<geometry::Sphere const&>(
+ 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()) {
+ auto const [t1, t2] = *opt;
+ if (t1.magnitude() >= 0)
+ intersectionTimes.push_back(t1);
+ else if (t2.magnitude() >= 0)
+ intersectionTimes.push_back(t2);
+ }
+ };
+
+ for (auto const& child : children) { addIfIntersects(*child); }
+
+ for (auto const* child : excluded) { addIfIntersects(*child); }
+
+ addIfIntersects(*currentVolumeNode);
+
+ auto const minIter =
+ std::min_element(intersectionTimes.cbegin(), intersectionTimes.cend());
+
+ TimeType min;
+
+ if (minIter == intersectionTimes.cend()) {
+ 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, min);
}
};
diff --git a/Processes/TrackingLine/testTrackingLine.cc b/Processes/TrackingLine/testTrackingLine.cc
new file mode 100644
index 0000000000000000000000000000000000000000..a3df326480c9417b736f6ac058dd70f6e8cce253
--- /dev/null
+++ b/Processes/TrackingLine/testTrackingLine.cc
@@ -0,0 +1,112 @@
+/**
+ * (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.
+ */
+
+#include <corsika/environment/Environment.h>
+
+#include <corsika/process/tracking_line/TrackingLine.h>
+
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/Sphere.h>
+#include <corsika/geometry/Vector.h>
+
+#include <corsika/setup/SetupTrajectory.h>
+using corsika::setup::Trajectory;
+
+#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
+ // cpp file
+#include <catch2/catch.hpp>
+
+using namespace corsika;
+using namespace corsika::process;
+using namespace corsika::units;
+using namespace corsika::geometry;
+
+#include <iostream>
+using namespace std;
+using namespace corsika::units::si;
+
+struct DummyParticle {
+ EnergyType fEnergy;
+ Vector<momentum_d> fMomentum;
+ Point fPosition;
+
+ DummyParticle(EnergyType pEnergy, Vector<momentum_d> pMomentum, Point pPosition)
+ : fEnergy(pEnergy)
+ , fMomentum(pMomentum)
+ , fPosition(pPosition) {}
+
+ auto GetEnergy() const { return fEnergy; }
+ auto GetMomentum() const { return fMomentum; }
+ auto GetPosition() const { return fPosition; }
+};
+
+struct DummyStack {
+ using ParticleType = DummyParticle;
+};
+
+TEST_CASE("TrackingLine") {
+ corsika::environment::Environment env; // dummy environment
+ auto const& cs = env.GetCoordinateSystem();
+
+ tracking_line::TrackingLine<DummyStack> tracking(env);
+
+ SECTION("intersection with sphere") {
+ Point const origin(cs, {0_m, 0_m, 0_m});
+ Point const center(cs, {0_m, 0_m, 10_m});
+ Sphere const sphere(center, 1_m);
+ Vector<corsika::units::si::SpeedType::dimension_type> v(cs, 0_m / second,
+ 0_m / second, 1_m / second);
+ Line line(origin, v);
+
+ geometry::Trajectory<Line> traj(line, 12345_s);
+
+ auto const opt =
+ tracking.TimeOfIntersection(traj, Sphere(Point(cs, {0_m, 0_m, 10_m}), 1_m));
+ REQUIRE(opt.has_value());
+
+ auto [t1, t2] = opt.value();
+ REQUIRE(t1 / 9_s == Approx(1));
+ REQUIRE(t2 / 11_s == Approx(1));
+
+ auto const optNoIntersection =
+ tracking.TimeOfIntersection(traj, Sphere(Point(cs, {5_m, 0_m, 10_m}), 1_m));
+ REQUIRE_FALSE(optNoIntersection.has_value());
+ }
+
+ SECTION("maximally possible propagation") {
+ auto& universe = *(env.GetUniverse());
+
+ //~ std::cout << env.GetUniverse().get() << std::endl;
+
+ DummyParticle p(
+ 1_J,
+ Vector<momentum_d>(cs, 0 * kilogram * meter / second,
+ 0 * kilogram * meter / second, 1 * kilogram * meter / second),
+ Point(cs, 0_m, 0_m, 0_m));
+
+ auto const radius = 20_m;
+
+ auto theMedium = corsika::environment::Environment::CreateNode<Sphere>(
+ Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m}, radius);
+ universe.AddChild(std::move(theMedium));
+
+ Point const origin(cs, {0_m, 0_m, 0_m});
+ Vector<corsika::units::si::SpeedType::dimension_type> v(cs, 0_m / second,
+ 0_m / second, 1_m / second);
+ Line line(origin, v);
+
+ auto const traj = tracking.GetTrack(p);
+
+ REQUIRE((traj.GetPosition(1.) - Point(cs, 0_m, 0_m, radius))
+ .GetComponents(cs)
+ .norm()
+ .magnitude() == Approx(0).margin(1e-4));
+ }
+}
diff --git a/Setup/SetupEnvironment.h b/Setup/SetupEnvironment.h
index 1beb2748052770fc25b84991fdf3e651744c3dc8..c5de9d53c9e26d6f5afe88d9762ec88df3043a64 100644
--- a/Setup/SetupEnvironment.h
+++ b/Setup/SetupEnvironment.h
@@ -12,6 +12,10 @@
#ifndef _include_corsika_setup_environment_h_
#define _include_corsika_setup_environment_h_
-namespace corsika {}
+#include <corsika/environment/IMediumModel.h>
+
+namespace corsika::setup {
+ using IEnvironmentModel = corsika::environment::IMediumModel;
+}
#endif
diff --git a/Stack/SuperStupidStack/SuperStupidStack.h b/Stack/SuperStupidStack/SuperStupidStack.h
index a1468d239e7f00227e79825b49439cdeba054e50..a889d32a652753955ba2553012350e521b1d99dc 100644
--- a/Stack/SuperStupidStack/SuperStupidStack.h
+++ b/Stack/SuperStupidStack/SuperStupidStack.h
@@ -128,7 +128,7 @@ namespace corsika::stack {
void Swap(const int i1, const int i2) {
std::swap(fDataPID[i2], fDataPID[i1]);
std::swap(fDataE[i2], fDataE[i1]);
- std::swap(fMomentum[i2], fMomentum[i1]); // should be Momentum !!!!
+ std::swap(fMomentum[i2], fMomentum[i1]);
std::swap(fPosition[i2], fPosition[i1]);
std::swap(fTime[i2], fTime[i1]);
}
@@ -137,7 +137,7 @@ namespace corsika::stack {
void IncrementSize() {
fDataPID.push_back(Code::Unknown);
fDataE.push_back(0 * joule);
- //#TODO this here makes no sense: see issue #48
+ //#TODO this here makes no sense: see issue #48
geometry::CoordinateSystem& dummyCS =
geometry::RootCoordinateSystem::GetInstance().GetRootCS();
fMomentum.push_back(MomentumVector(