diff --git a/Documentation/Examples/CMakeLists.txt b/Documentation/Examples/CMakeLists.txt
index f9c631669868f15f8d2d1ea55dcf18ed970a96e1..c677f3ab9a3dccd82235fee8d4690be4de697d17 100644
--- a/Documentation/Examples/CMakeLists.txt
+++ b/Documentation/Examples/CMakeLists.txt
@@ -22,30 +22,27 @@ target_link_libraries (stack_example SuperStupidStack CORSIKAunits
CORSIKAlogging)
CORSIKA_ADD_TEST (stack_example)
-if (Pythia8_FOUND)
- add_executable (cascade_example cascade_example.cc)
- target_compile_options(cascade_example PRIVATE -g) # do not skip asserts
- target_link_libraries (cascade_example SuperStupidStack CORSIKAunits
- CORSIKAlogging
- CORSIKArandom
- ProcessSibyll
- ProcessPythia
- CORSIKAcascade
- ProcessParticleCut
- ProcessEnergyLoss
- ProcessStackInspector
- ProcessTrackWriter
- ProcessTrackingLine
- CORSIKAprocesses
- CORSIKAcascade
- CORSIKAparticles
- CORSIKAgeometry
- CORSIKAenvironment
- CORSIKAprocesssequence
- )
- install (TARGETS cascade_example DESTINATION share/examples)
- CORSIKA_ADD_TEST (cascade_example)
-endif()
+add_executable (cascade_example cascade_example.cc)
+target_compile_options(cascade_example PRIVATE -g) # do not skip asserts
+target_link_libraries (cascade_example SuperStupidStack CORSIKAunits
+ CORSIKAlogging
+ CORSIKArandom
+ ProcessSibyll
+ CORSIKAcascade
+ ProcessEnergyLoss
+ ProcessStackInspector
+ ProcessParticleCut
+ ProcessTrackWriter
+ ProcessTrackingLine
+ CORSIKAprocesses
+ CORSIKAcascade
+ CORSIKAparticles
+ CORSIKAgeometry
+ CORSIKAenvironment
+ CORSIKAprocesssequence
+ )
+install (TARGETS cascade_example DESTINATION share/examples)
+CORSIKA_ADD_TEST (cascade_example)
add_executable (boundary_example boundary_example.cc)
target_compile_options(boundary_example PRIVATE -g) # do not skip asserts
@@ -62,7 +59,6 @@ target_link_libraries (boundary_example SuperStupidStack CORSIKAunits CORSIKAlog
CORSIKAenvironment
CORSIKAprocesssequence
)
-
install (TARGETS boundary_example DESTINATION share/examples)
CORSIKA_ADD_TEST (boundary_example)
@@ -80,6 +76,7 @@ if (Pythia8_FOUND)
ProcessTrackingLine
ProcessParticleCut
ProcessHadronicElasticModel
+ ProcessStackInspector
CORSIKAprocesses
CORSIKAcascade
CORSIKAparticles
@@ -92,30 +89,28 @@ if (Pythia8_FOUND)
CORSIKA_ADD_TEST (cascade_proton_example)
endif()
-if (Pythia8_FOUND)
- add_executable (vertical_EAS vertical_EAS.cc)
- target_compile_options(vertical_EAS PRIVATE -g) # do not skip asserts
- target_link_libraries (vertical_EAS SuperStupidStack CORSIKAunits
- CORSIKAlogging
- CORSIKArandom
- ProcessSibyll
- ProcessPythia
- CORSIKAcascade
- ProcessEnergyLoss
- ProcessTrackWriter
- ProcessTrackingLine
- ProcessParticleCut
- ProcessHadronicElasticModel
- CORSIKAprocesses
- CORSIKAcascade
- CORSIKAparticles
- CORSIKAgeometry
- CORSIKAenvironment
- CORSIKAprocesssequence
- )
- install (TARGETS vertical_EAS DESTINATION share/examples)
- # CORSIKA_ADD_TEST (vertical_EAS) not yet...
-endif()
+add_executable (vertical_EAS vertical_EAS.cc)
+target_compile_options(vertical_EAS PRIVATE -g) # do not skip asserts
+target_link_libraries (vertical_EAS SuperStupidStack CORSIKAunits
+ CORSIKAlogging
+ CORSIKArandom
+ ProcessSibyll
+ CORSIKAcascade
+ ProcessEnergyLoss
+ ProcessTrackWriter
+ ProcessTrackingLine
+ ProcessHadronicElasticModel
+ ProcessParticleCut
+ ProcessStackInspector
+ CORSIKAprocesses
+ CORSIKAcascade
+ CORSIKAparticles
+ CORSIKAgeometry
+ CORSIKAenvironment
+ CORSIKAprocesssequence
+ )
+install (TARGETS vertical_EAS DESTINATION share/examples)
+# CORSIKA_ADD_TEST (vertical_EAS) not yet...
add_executable (staticsequence_example staticsequence_example.cc)
target_compile_options(staticsequence_example PRIVATE -g) # do not skip asserts
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 400bb07d8166075b038cd8b59ded4ae759769808..7690d664811f4192d224044ab39ebea7b92d6136 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -29,8 +29,6 @@
#include <corsika/process/sibyll/Interaction.h>
#include <corsika/process/sibyll/NuclearInteraction.h>
-#include <corsika/process/pythia/Decay.h>
-
#include <corsika/process/track_writer/TrackWriter.h>
#include <corsika/process/particle_cut/ParticleCut.h>
@@ -105,7 +103,6 @@ int main() {
particles::Code::KMinus, particles::Code::K0Long, particles::Code::K0Short};
random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
- random::RNGManager::GetInstance().RegisterRandomStream("pythia");
process::sibyll::Interaction sibyll;
process::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
process::sibyll::Decay decay(trackedHadrons);
diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc
index f96544ffbc5c8034f3e06b711590e2885f90af4c..872ea2d943bb808edeb295b767cc931dc3e49032 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -28,10 +28,7 @@
#include <corsika/process/sibyll/Interaction.h>
#include <corsika/process/sibyll/NuclearInteraction.h>
-#include <corsika/process/pythia/Decay.h>
-
#include <corsika/process/particle_cut/ParticleCut.h>
-
#include <corsika/process/track_writer/TrackWriter.h>
#include <corsika/units/PhysicalUnits.h>
@@ -59,7 +56,6 @@ using namespace corsika::environment;
using namespace std;
using namespace corsika::units::si;
-
//
// The example main program for a particle cascade
//
@@ -102,8 +98,6 @@ int main() {
process::sibyll::Interaction sibyll;
process::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
process::sibyll::Decay decay(trackedHadrons);
- // random::RNGManager::GetInstance().RegisterRandomStream("pythia");
- // process::pythia::Decay decay(trackedHadrons);
process::particle_cut::ParticleCut cut(20_GeV);
// random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel");
// process::HadronicElasticModel::HadronicElasticInteraction
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index 75f29e4ed51bc70d7d31d4917ef4c186baee4bef..f9ec0beb0e40df08a16595bee6be8dc0048cff94 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -153,7 +153,7 @@ namespace corsika::cascade {
// determine combined total interaction length (inverse)
InverseGrammageType const total_inv_lambda =
- fProcessSequence.GetTotalInverseInteractionLength(vParticle, step);
+ fProcessSequence.GetTotalInverseInteractionLength(vParticle);
// sample random exponential step length in grammage
corsika::random::ExponentialDistribution expDist(1 / total_inv_lambda);
@@ -242,13 +242,14 @@ namespace corsika::cascade {
if (min_distance == distance_interact) {
std::cout << "collide" << std::endl;
- InverseGrammageType const actual_inv_length =
- fProcessSequence.GetTotalInverseInteractionLength(vParticle, step);
+ InverseGrammageType const current_inv_length =
+ fProcessSequence.GetTotalInverseInteractionLength(vParticle);
- random::UniformRealDistribution<InverseGrammageType> uniDist(actual_inv_length);
+ random::UniformRealDistribution<InverseGrammageType> uniDist(
+ current_inv_length);
const auto sample_process = uniDist(fRNG);
InverseGrammageType inv_lambda_count = 0. * meter * meter / gram;
- fProcessSequence.SelectInteraction(vParticle, projectile, step, sample_process,
+ fProcessSequence.SelectInteraction(vParticle, projectile, sample_process,
inv_lambda_count);
} else if (min_distance == distance_decay) {
std::cout << "decay" << std::endl;
diff --git a/Framework/Cascade/testCascade.cc b/Framework/Cascade/testCascade.cc
index 89fd65ef42150149a9c0ffbd7183fed865ba659d..f29d1af9fba4c01c6452fc2dc026df7ff99f2c42 100644
--- a/Framework/Cascade/testCascade.cc
+++ b/Framework/Cascade/testCascade.cc
@@ -69,8 +69,8 @@ public:
ProcessSplit(GrammageType const X0)
: fX0(X0) {}
- template <typename Particle, typename Track>
- corsika::units::si::GrammageType GetInteractionLength(Particle&, Track&) const {
+ template <typename Particle>
+ corsika::units::si::GrammageType GetInteractionLength(Particle const&) const {
return fX0;
}
diff --git a/Framework/Geometry/CMakeLists.txt b/Framework/Geometry/CMakeLists.txt
index c12c79b01d9663c8caf08a8068ec0dc13b529e3c..510d8e3922f2c7c9a1bbd0578bea34f74641e500 100644
--- a/Framework/Geometry/CMakeLists.txt
+++ b/Framework/Geometry/CMakeLists.txt
@@ -10,6 +10,7 @@ set (
Point.h
Line.h
Sphere.h
+ Plane.h
Volume.h
CoordinateSystem.h
RootCoordinateSystem.h
diff --git a/Framework/Geometry/Plane.h b/Framework/Geometry/Plane.h
new file mode 100644
index 0000000000000000000000000000000000000000..2c3f103f52bb9278ffcca63766512ceda13d32c2
--- /dev/null
+++ b/Framework/Geometry/Plane.h
@@ -0,0 +1,41 @@
+
+/*
+ * (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_Framework_Geometry_Plane_h_
+#define _include_Framework_Geometry_Plane_h_
+
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/Vector.h>
+#include <corsika/units/PhysicalUnits.h>
+
+namespace corsika::geometry {
+ class Plane {
+
+ using DimLessVec = Vector<corsika::units::si::dimensionless_d>;
+
+ Point const fCenter;
+ DimLessVec const fNormal;
+
+ public:
+ Plane(Point const& vCenter, DimLessVec const& vNormal)
+ : fCenter(vCenter)
+ , fNormal(vNormal) {}
+ bool IsAbove(Point const& vP) const {
+ return fNormal.dot(vP - fCenter) > corsika::units::si::LengthType::zero();
+ }
+
+ Point const& GetCenter() const { return fCenter; }
+ DimLessVec const& GetNormal() const { return fNormal; }
+ };
+
+} // namespace corsika::geometry
+
+#endif
diff --git a/Framework/ProcessSequence/ContinuousProcess.h b/Framework/ProcessSequence/ContinuousProcess.h
index 58d378aa3eb5bdad5f8029555ea8825afa450007..27f6c8fae4cb52e037f7ac128dd5fde161dd5566 100644
--- a/Framework/ProcessSequence/ContinuousProcess.h
+++ b/Framework/ProcessSequence/ContinuousProcess.h
@@ -34,11 +34,11 @@ namespace corsika::process {
// here starts the interface part
// -> enforce derived to implement DoContinuous...
template <typename Particle, typename Track>
- EProcessReturn DoContinuous(Particle&, Track&) const;
+ EProcessReturn DoContinuous(Particle&, Track const&) const;
// -> enforce derived to implement MaxStepLength...
template <typename Particle, typename Track>
- corsika::units::si::LengthType MaxStepLength(Particle& p, Track& track) const;
+ units::si::LengthType MaxStepLength(Particle const& p, Track const& track) const;
};
// overwrite the default trait class, to mark BaseProcess<T> as useful process
diff --git a/Framework/ProcessSequence/InteractionProcess.h b/Framework/ProcessSequence/InteractionProcess.h
index 27ad9c222d5ce644d6b344a469edaf8e717fcf87..a9b4b9fca306e865a85ac94ae6a9878ce7f0057d 100644
--- a/Framework/ProcessSequence/InteractionProcess.h
+++ b/Framework/ProcessSequence/InteractionProcess.h
@@ -38,13 +38,12 @@ namespace corsika::process {
template <typename Particle>
EProcessReturn DoInteraction(Particle&);
- template <typename Particle, typename Track>
- corsika::units::si::GrammageType GetInteractionLength(Particle& p, Track& t);
+ template <typename TParticle>
+ corsika::units::si::GrammageType GetInteractionLength(TParticle& p);
- template <typename Particle, typename Track>
- corsika::units::si::InverseGrammageType GetInverseInteractionLength(Particle& p,
- Track& t) {
- return 1. / GetRef().GetInteractionLength(p, t);
+ template <typename TParticle>
+ corsika::units::si::InverseGrammageType GetInverseInteractionLength(TParticle& p) {
+ return 1. / GetRef().GetInteractionLength(p);
}
};
diff --git a/Framework/ProcessSequence/ProcessSequence.h b/Framework/ProcessSequence/ProcessSequence.h
index 3bd759c04a97e55360c87473cc4457159bd967fa..55eb43ba795c35d04654a5256232221fa4c900b5 100644
--- a/Framework/ProcessSequence/ProcessSequence.h
+++ b/Framework/ProcessSequence/ProcessSequence.h
@@ -47,9 +47,19 @@ namespace corsika::process {
// this is a marker to track which BaseProcess is also a ProcessSequence
template <typename T>
- struct is_process_sequence {
- static const bool value = false;
- };
+ struct is_process_sequence : std::false_type {};
+
+ template <typename T>
+ bool constexpr is_process_sequence_v = is_process_sequence<T>::value;
+
+ namespace switch_process {
+ template <typename A, typename B>
+ class SwitchProcess; // fwd-decl.
+ }
+
+ template <typename A, typename B>
+ struct is_process_sequence<switch_process::SwitchProcess<A, B>>
+ : std::true_type {};
/**
T1 and T2 are both references if possible (lvalue), otherwise
@@ -57,11 +67,14 @@ namespace corsika::process {
rvalue as well as lvalue Processes in the ProcessSequence.
*/
template <typename T1, typename T2>
- class ProcessSequence : public BaseProcess<ProcessSequence<T1, T2> > {
+ class ProcessSequence : public BaseProcess<ProcessSequence<T1, T2>> {
using T1type = typename std::decay<T1>::type;
using T2type = typename std::decay<T2>::type;
+ static bool constexpr t1ProcSeq = is_process_sequence_v<T1type>;
+ static bool constexpr t2ProcSeq = is_process_sequence_v<T2type>;
+
public:
T1 A; // this is a reference, if possible
T2 B; // this is a reference, if possible
@@ -78,13 +91,13 @@ namespace corsika::process {
VTNType const& to) {
EProcessReturn ret = EProcessReturn::eOk;
- if constexpr (std::is_base_of<BoundaryCrossingProcess<T1type>, T1type>::value ||
- is_process_sequence<T1>::value) {
+ if constexpr (std::is_base_of_v<BoundaryCrossingProcess<T1type>, T1type> ||
+ t1ProcSeq) {
ret |= A.DoBoundaryCrossing(p, from, to);
}
- if constexpr (std::is_base_of<BoundaryCrossingProcess<T2type>, T2type>::value ||
- is_process_sequence<T2>::value) {
+ if constexpr (std::is_base_of_v<BoundaryCrossingProcess<T2type>, T2type> ||
+ t2ProcSeq) {
ret |= B.DoBoundaryCrossing(p, from, to);
}
@@ -94,12 +107,10 @@ namespace corsika::process {
template <typename TParticle, typename TTrack>
EProcessReturn DoContinuous(TParticle& vP, TTrack& vT) {
EProcessReturn ret = EProcessReturn::eOk;
- if constexpr (std::is_base_of<ContinuousProcess<T1type>, T1type>::value ||
- is_process_sequence<T1>::value) {
+ if constexpr (std::is_base_of_v<ContinuousProcess<T1type>, T1type> || t1ProcSeq) {
ret |= A.DoContinuous(vP, vT);
}
- if constexpr (std::is_base_of<ContinuousProcess<T2type>, T2type>::value ||
- is_process_sequence<T2>::value) {
+ if constexpr (std::is_base_of_v<ContinuousProcess<T2type>, T2type> || t2ProcSeq) {
ret |= B.DoContinuous(vP, vT);
}
return ret;
@@ -108,12 +119,10 @@ namespace corsika::process {
template <typename TSecondaries>
EProcessReturn DoSecondaries(TSecondaries& vS) {
EProcessReturn ret = EProcessReturn::eOk;
- if constexpr (std::is_base_of<SecondariesProcess<T1type>, T1type>::value ||
- is_process_sequence<T1>::value) {
+ if constexpr (std::is_base_of_v<SecondariesProcess<T1type>, T1type> || t1ProcSeq) {
ret |= A.DoSecondaries(vS);
}
- if constexpr (std::is_base_of<SecondariesProcess<T2type>, T2type>::value ||
- is_process_sequence<T2>::value) {
+ if constexpr (std::is_base_of_v<SecondariesProcess<T2type>, T2type> || t2ProcSeq) {
ret |= B.DoSecondaries(vS);
}
return ret;
@@ -121,12 +130,10 @@ namespace corsika::process {
bool CheckStep() {
bool ret = false;
- if constexpr (std::is_base_of<StackProcess<T1type>, T1type>::value ||
- is_process_sequence<T1>::value) {
+ if constexpr (std::is_base_of_v<StackProcess<T1type>, T1type> || t1ProcSeq) {
ret |= A.CheckStep();
}
- if constexpr (std::is_base_of<StackProcess<T2type>, T2type>::value ||
- is_process_sequence<T2>::value) {
+ if constexpr (std::is_base_of_v<StackProcess<T2type>, T2type> || t2ProcSeq) {
ret |= B.CheckStep();
}
return ret;
@@ -135,12 +142,10 @@ namespace corsika::process {
template <typename TStack>
EProcessReturn DoStack(TStack& vS) {
EProcessReturn ret = EProcessReturn::eOk;
- if constexpr (std::is_base_of<StackProcess<T1type>, T1type>::value ||
- is_process_sequence<T1>::value) {
+ if constexpr (std::is_base_of_v<StackProcess<T1type>, T1type> || t1ProcSeq) {
if (A.CheckStep()) { ret |= A.DoStack(vS); }
}
- if constexpr (std::is_base_of<StackProcess<T2type>, T2type>::value ||
- is_process_sequence<T2>::value) {
+ if constexpr (std::is_base_of_v<StackProcess<T2type>, T2type> || t2ProcSeq) {
if (B.CheckStep()) { ret |= B.DoStack(vS); }
}
return ret;
@@ -152,64 +157,58 @@ namespace corsika::process {
max_length = // if no other process in the sequence implements it
std::numeric_limits<double>::infinity() * corsika::units::si::meter;
- if constexpr (std::is_base_of<ContinuousProcess<T1type>, T1type>::value ||
- is_process_sequence<T1>::value) {
+ if constexpr (std::is_base_of_v<ContinuousProcess<T1type>, T1type> || t1ProcSeq) {
corsika::units::si::LengthType const len = A.MaxStepLength(vP, vTrack);
max_length = std::min(max_length, len);
}
- if constexpr (std::is_base_of<ContinuousProcess<T2type>, T2type>::value ||
- is_process_sequence<T2>::value) {
+ if constexpr (std::is_base_of_v<ContinuousProcess<T2type>, T2type> || t2ProcSeq) {
corsika::units::si::LengthType const len = B.MaxStepLength(vP, vTrack);
max_length = std::min(max_length, len);
}
return max_length;
}
- template <typename TParticle, typename TTrack>
- corsika::units::si::GrammageType GetTotalInteractionLength(TParticle& vP,
- TTrack& vT) {
- return 1. / GetInverseInteractionLength(vP, vT);
+ template <typename TParticle>
+ corsika::units::si::GrammageType GetTotalInteractionLength(TParticle& vP) {
+ return 1. / GetInverseInteractionLength(vP);
}
- template <typename TParticle, typename TTrack>
+ template <typename TParticle>
corsika::units::si::InverseGrammageType GetTotalInverseInteractionLength(
- TParticle& vP, TTrack& vT) {
- return GetInverseInteractionLength(vP, vT);
+ TParticle& vP) {
+ return GetInverseInteractionLength(vP);
}
- template <typename TParticle, typename TTrack>
- corsika::units::si::InverseGrammageType GetInverseInteractionLength(TParticle& vP,
- TTrack& vT) {
+ template <typename TParticle>
+ corsika::units::si::InverseGrammageType GetInverseInteractionLength(TParticle& vP) {
using namespace corsika::units::si;
InverseGrammageType tot = 0 * meter * meter / gram;
- if constexpr (std::is_base_of<InteractionProcess<T1type>, T1type>::value ||
- is_process_sequence<T1>::value) {
- tot += A.GetInverseInteractionLength(vP, vT);
+ if constexpr (std::is_base_of_v<InteractionProcess<T1type>, T1type> || t1ProcSeq) {
+ tot += A.GetInverseInteractionLength(vP);
}
- if constexpr (std::is_base_of<InteractionProcess<T2type>, T2type>::value ||
- is_process_sequence<T2>::value) {
- tot += B.GetInverseInteractionLength(vP, vT);
+ if constexpr (std::is_base_of_v<InteractionProcess<T2type>, T2type> || t2ProcSeq) {
+ tot += B.GetInverseInteractionLength(vP);
}
return tot;
}
- template <typename TParticle, typename TSecondaries, typename TTrack>
+ template <typename TParticle, typename TSecondaries>
EProcessReturn SelectInteraction(
- TParticle& vP, TSecondaries& vS, TTrack& vT,
+ TParticle& vP, TSecondaries& vS,
[[maybe_unused]] corsika::units::si::InverseGrammageType lambda_select,
corsika::units::si::InverseGrammageType& lambda_inv_count) {
- if constexpr (is_process_sequence<T1type>::value) {
+ if constexpr (t1ProcSeq) {
// if A is a process sequence --> check inside
const EProcessReturn ret =
- A.SelectInteraction(vP, vS, vT, lambda_select, lambda_inv_count);
+ A.SelectInteraction(vP, vS, lambda_select, lambda_inv_count);
// if A did succeed, stop routine
if (ret != EProcessReturn::eOk) { return ret; }
- } else if constexpr (std::is_base_of<InteractionProcess<T1type>, T1type>::value) {
+ } else if constexpr (std::is_base_of_v<InteractionProcess<T1type>, T1type>) {
// if this is not a ContinuousProcess --> evaluate probability
- lambda_inv_count += A.GetInverseInteractionLength(vP, vT);
+ lambda_inv_count += A.GetInverseInteractionLength(vP);
// check if we should execute THIS process and then EXIT
if (lambda_select < lambda_inv_count) {
A.DoInteraction(vS);
@@ -217,15 +216,15 @@ namespace corsika::process {
}
} // end branch A
- if constexpr (is_process_sequence<T2>::value) {
+ if constexpr (t2ProcSeq) {
// if A is a process sequence --> check inside
const EProcessReturn ret =
- B.SelectInteraction(vP, vS, vT, lambda_select, lambda_inv_count);
+ B.SelectInteraction(vP, vS, lambda_select, lambda_inv_count);
// if A did succeed, stop routine
if (ret != EProcessReturn::eOk) { return ret; }
- } else if constexpr (std::is_base_of<InteractionProcess<T2type>, T2type>::value) {
+ } else if constexpr (std::is_base_of_v<InteractionProcess<T2type>, T2type>) {
// if this is not a ContinuousProcess --> evaluate probability
- lambda_inv_count += B.GetInverseInteractionLength(vP, vT);
+ lambda_inv_count += B.GetInverseInteractionLength(vP);
// check if we should execute THIS process and then EXIT
if (lambda_select < lambda_inv_count) {
B.DoInteraction(vS);
@@ -251,12 +250,10 @@ namespace corsika::process {
corsika::units::si::InverseTimeType tot = 0 / second;
- if constexpr (std::is_base_of<DecayProcess<T1type>, T1type>::value ||
- is_process_sequence<T1>::value) {
+ if constexpr (std::is_base_of_v<DecayProcess<T1type>, T1type> || t1ProcSeq) {
tot += A.GetInverseLifetime(p);
}
- if constexpr (std::is_base_of<DecayProcess<T2type>, T2type>::value ||
- is_process_sequence<T2>::value) {
+ if constexpr (std::is_base_of_v<DecayProcess<T2type>, T2type> || t2ProcSeq) {
tot += B.GetInverseLifetime(p);
}
return tot;
@@ -268,12 +265,12 @@ namespace corsika::process {
TParticle& vP, TSecondaries& vS,
[[maybe_unused]] corsika::units::si::InverseTimeType decay_select,
corsika::units::si::InverseTimeType& decay_inv_count) {
- if constexpr (is_process_sequence<T1>::value) {
+ if constexpr (t1ProcSeq) {
// if A is a process sequence --> check inside
const EProcessReturn ret = A.SelectDecay(vP, vS, decay_select, decay_inv_count);
// if A did succeed, stop routine
if (ret != EProcessReturn::eOk) { return ret; }
- } else if constexpr (std::is_base_of<DecayProcess<T1type>, T1type>::value) {
+ } else if constexpr (std::is_base_of_v<DecayProcess<T1type>, T1type>) {
// if this is not a ContinuousProcess --> evaluate probability
decay_inv_count += A.GetInverseLifetime(vP);
// check if we should execute THIS process and then EXIT
@@ -284,12 +281,12 @@ namespace corsika::process {
}
} // end branch A
- if constexpr (is_process_sequence<T2>::value) {
+ if constexpr (t2ProcSeq) {
// if A is a process sequence --> check inside
const EProcessReturn ret = B.SelectDecay(vP, vS, decay_select, decay_inv_count);
// if A did succeed, stop routine
if (ret != EProcessReturn::eOk) { return ret; }
- } else if constexpr (std::is_base_of<DecayProcess<T2type>, T2type>::value) {
+ } else if constexpr (std::is_base_of_v<DecayProcess<T2type>, T2type>) {
// if this is not a ContinuousProcess --> evaluate probability
decay_inv_count += B.GetInverseLifetime(vP);
// check if we should execute THIS process and then EXIT
@@ -322,9 +319,7 @@ namespace corsika::process {
/// marker to identify objectas ProcessSequence
template <typename A, typename B>
- struct is_process_sequence<corsika::process::ProcessSequence<A, B> > {
- static const bool value = true;
- };
+ struct is_process_sequence<corsika::process::ProcessSequence<A, B>> : std::true_type {};
} // namespace corsika::process
diff --git a/Framework/ProcessSequence/testProcessSequence.cc b/Framework/ProcessSequence/testProcessSequence.cc
index aa78ee89c4c012eb2df3d322b6b1ddf6fa44ba1e..e8bfe7c1c835762e64ccb3afed48962775dc24ae 100644
--- a/Framework/ProcessSequence/testProcessSequence.cc
+++ b/Framework/ProcessSequence/testProcessSequence.cc
@@ -100,8 +100,8 @@ public:
cout << "Process2::DoInteraction" << endl;
return EProcessReturn::eOk;
}
- template <typename Particle, typename Track>
- GrammageType GetInteractionLength(Particle&, Track&) const {
+ template <typename Particle>
+ GrammageType GetInteractionLength(Particle&) const {
cout << "Process2::GetInteractionLength" << endl;
return 3_g / (1_cm * 1_cm);
}
@@ -123,8 +123,8 @@ public:
cout << "Process3::DoInteraction" << endl;
return EProcessReturn::eOk;
}
- template <typename Particle, typename Track>
- GrammageType GetInteractionLength(Particle&, Track&) const {
+ template <typename Particle>
+ GrammageType GetInteractionLength(Particle&) const {
cout << "Process3::GetInteractionLength" << endl;
return 1_g / (1_cm * 1_cm);
}
@@ -220,12 +220,11 @@ TEST_CASE("Process Sequence", "[Process Sequence]") {
Process3 m3(2);
DummyData particle;
- DummyTrajectory track;
auto sequence2 = cp1 << m2 << m3;
- GrammageType const tot = sequence2.GetTotalInteractionLength(particle, track);
+ GrammageType const tot = sequence2.GetTotalInteractionLength(particle);
InverseGrammageType const tot_inv =
- sequence2.GetTotalInverseInteractionLength(particle, track);
+ sequence2.GetTotalInverseInteractionLength(particle);
cout << "lambda_tot=" << tot << "; lambda_tot_inv=" << tot_inv << endl;
}
diff --git a/Framework/StackInterface/SecondaryView.h b/Framework/StackInterface/SecondaryView.h
index 1bbbe02ccaad07fcc533d5782c753427c0dddbf1..fae9d463f9c6808b69d346092fe6277a7c5fec7a 100644
--- a/Framework/StackInterface/SecondaryView.h
+++ b/Framework/StackInterface/SecondaryView.h
@@ -43,14 +43,18 @@ namespace corsika::stack {
All operations possible on a Stack object are also possible on a
SecondaryView object. This means you can add, delete, copy, swap,
iterate, etc.
- */
- /* INTERNAL NOTE FOR DEVELOPERS The secondary particle indices are
- stored in a std::vector fIndices, the index of the primary projectle
- particle is explicitly stored in fProjectileIndex. StackIterator
- indices are refering to those numbers, where
- StackIterator::GetIndex()==0 refers to the fProjectileIndex and
- StackIterator::GetIndex()>0 to fIndices[i+1], see function GetIndexFromIterator.
+ *Further information about implementation (for developers):* All
+ data is stored in the original stack privided at construction
+ time. The secondary particle (view) indices are stored in an
+ extra std::vector of SecondaryView class 'fIndices' referring to
+ the original stack slot indices. The index of the primary
+ projectle particle is also explicitly stored in
+ 'fProjectileIndex'. StackIterator indices
+ 'i = StackIterator::GetIndex()' are referring to those numbers,
+ where 'i==0' refers to the 'fProjectileIndex', and
+ 'StackIterator::GetIndex()>0' to 'fIndices[i-1]', see function
+ GetIndexFromIterator.
*/
template <typename StackDataType, template <typename> typename ParticleInterface>
@@ -125,8 +129,11 @@ namespace corsika::stack {
template <typename... Args>
auto AddSecondary(StackIterator& proj, const Args... v) {
+ // make space on stack
InnerStackType::GetStackData().IncrementSize();
+ // get current number of secondaries on stack
const unsigned int idSec = GetSize();
+ // determine index on (inner) stack where new particle will be located
const unsigned int index = InnerStackType::GetStackData().GetSize() - 1;
fIndices.push_back(index);
// NOTE: "+1" is since "0" is special marker here for PROJECTILE, see
@@ -161,14 +168,26 @@ namespace corsika::stack {
/// @}
/**
- * need overwrite Stack::Delete, since we want to call SecondaryView::DeleteLast
+ * need overwrite Stack::Delete, since we want to call
+ * SecondaryView::DeleteLast
+ *
+ * The particle is deleted on the underlying (internal) stack. The
+ * local references in SecondaryView in fIndices must be fixed,
+ * too. The approach is to a) check if the particle 'p' is at the
+ * very end of the internal stack, b) if not: move it there by
+ * copying the last particle to the current particle location, c)
+ * remove the last particle.
+ *
*/
void Delete(StackIterator p) {
if (IsEmpty()) { /* error */
throw std::runtime_error("Stack, cannot delete entry since size is zero");
}
- if (p.GetIndex() < GetSize() - 1)
- InnerStackType::GetStackData().Copy(GetSize() - 1, p.GetIndex());
+ const int innerSize = InnerStackType::GetSize();
+ const int innerIndex = GetIndexFromIterator(p.GetIndex());
+ if (innerIndex < innerSize - 1)
+ InnerStackType::GetStackData().Copy(innerSize - 1,
+ GetIndexFromIterator(p.GetIndex()));
DeleteLast();
}
diff --git a/Framework/StackInterface/testCombinedStack.cc b/Framework/StackInterface/testCombinedStack.cc
index 8fe244b667bf736cc418368ee2e06116ef845608..aceccd1b472a99c57168eef4d78851b3641c092c 100644
--- a/Framework/StackInterface/testCombinedStack.cc
+++ b/Framework/StackInterface/testCombinedStack.cc
@@ -333,7 +333,8 @@ using StackTest2 = CombinedStack<typename StackTest::StackImpl, TestStackData3,
CombinedTestInterfaceType2>;
#if defined(__clang__)
-using StackTestView = SecondaryView<typename StackTest2::StackImpl, CombinedTestInterfaceType2>;
+using StackTestView =
+ SecondaryView<typename StackTest2::StackImpl, CombinedTestInterfaceType2>;
#elif defined(__GNUC__) || defined(__GNUG__)
using StackTestView = corsika::stack::MakeView<StackTest2>::type;
#endif
diff --git a/Framework/StackInterface/testSecondaryView.cc b/Framework/StackInterface/testSecondaryView.cc
index 2d2064b64e2e09c96ec5f77ebbf25d50f5d60e9e..95c196d5ecef139339b21bb41623f8a52e846d76 100644
--- a/Framework/StackInterface/testSecondaryView.cc
+++ b/Framework/StackInterface/testSecondaryView.cc
@@ -134,4 +134,67 @@ TEST_CASE("SecondaryStack", "[stack]") {
REQUIRE(view.GetSize() == 1);
}
+
+ SECTION("deletion") {
+ StackTest stack;
+ stack.AddParticle(std::tuple{-99.});
+ stack.AddParticle(std::tuple{0.});
+
+ {
+ auto particle = stack.GetNextParticle();
+ StackTestView view(particle);
+
+ auto proj = view.GetProjectile();
+ proj.AddSecondary(std::tuple{-2.});
+ proj.AddSecondary(std::tuple{-1.});
+ proj.AddSecondary(std::tuple{1.});
+ proj.AddSecondary(std::tuple{2.});
+
+ CHECK(stack.GetSize() == 6); // -99, 0, -2, -1, 1, 2
+ CHECK(view.GetSize() == 4); // -2, -1, 1, 2
+
+ // now delete all negative entries, i.e. -1 and -2
+ auto p = view.begin();
+ while (p != view.end()) {
+ auto data = p.GetData();
+ if (data < 0) {
+ p.Delete();
+ } else {
+ ++p;
+ }
+ }
+ CHECK(stack.GetSize() == 4); // -99, 0, 2, 1 (order changes during deletion)
+ CHECK(view.GetSize() == 2); // 2, 1
+ }
+
+ // repeat
+
+ {
+ auto particle = stack.GetNextParticle();
+ StackTestView view(particle);
+
+ // put -2,...,+2 on stack
+ auto proj = view.GetProjectile();
+ proj.AddSecondary(std::tuple{-2.});
+ proj.AddSecondary(std::tuple{-1.});
+ proj.AddSecondary(std::tuple{1.});
+ proj.AddSecondary(std::tuple{2.});
+ // stack should contain -99, 0, 2, 1, [-2, -1, 1, 2]
+
+ auto p = view.begin();
+ while (p != view.end()) {
+ auto data = p.GetData();
+ if (data < 0) {
+ p.Delete();
+ } else {
+ ++p;
+ }
+ }
+
+ // stack should contain -99, 0, 2, 1, [2, 1]
+ // view should contain 1, 2
+
+ CHECK(stack.GetSize() == 6);
+ }
+ }
}
diff --git a/Framework/Units/PhysicalConstants.h b/Framework/Units/PhysicalConstants.h
index 75dfc51cedd83bd266f07f971aef48ca717b7d81..23494f036ed2f692de71984c8fc1b81a6d919112 100644
--- a/Framework/Units/PhysicalConstants.h
+++ b/Framework/Units/PhysicalConstants.h
@@ -51,6 +51,9 @@ namespace corsika::units::constants {
auto constexpr nucleonMass = 0.5 * (0.93827 + 0.93957) * 1e9 * electronvolt;
+ // molar gas constant
+ auto constexpr R = Rep(8.314'459'8) * joule / (mole * kelvin);
+
// etc.
} // namespace corsika::units::constants
diff --git a/Framework/Utilities/sgn.h b/Framework/Utilities/sgn.h
index 71d5bec5bbb59fe8c71ef35d0c3f9da1a40b27b7..43d9a9fad728785ec851424b2113f5e7dac60af1 100644
--- a/Framework/Utilities/sgn.h
+++ b/Framework/Utilities/sgn.h
@@ -3,12 +3,12 @@
namespace corsika::utl {
-//! sign function without branches
-template <typename T>
-static int sgn(T val) {
- return (T(0) < val) - (val < T(0));
-}
+ //! sign function without branches
+ template <typename T>
+ static int sgn(T val) {
+ return (T(0) < val) - (val < T(0));
+ }
-}
+} // namespace corsika::utl
#endif
diff --git a/Processes/CMakeLists.txt b/Processes/CMakeLists.txt
index c9f218bf0066f3b7d16f883464090beca71e8000..b54f76d8fd3d16d37de019ffb202015ab15bd663 100644
--- a/Processes/CMakeLists.txt
+++ b/Processes/CMakeLists.txt
@@ -10,6 +10,7 @@ add_subdirectory (TrackingLine)
add_subdirectory (HadronicElasticModel)
add_subdirectory (EnergyLoss)
add_subdirectory (UrQMD)
+add_subdirectory (SwitchProcess)
add_subdirectory (ParticleCut)
#add_custom_target(CORSIKAprocesses)
diff --git a/Processes/EnergyLoss/EnergyLoss.cc b/Processes/EnergyLoss/EnergyLoss.cc
index 990d6172501dc67f0df4958cf651d587003f1923..aaa6dfd6146deb0f34cfbdd79f7e78d76aafa4be 100644
--- a/Processes/EnergyLoss/EnergyLoss.cc
+++ b/Processes/EnergyLoss/EnergyLoss.cc
@@ -24,9 +24,8 @@ using namespace std;
using namespace corsika;
using namespace corsika::units::si;
-using namespace corsika::setup;
-using Particle = Stack::ParticleType;
-using Track = Trajectory;
+using SetupParticle = corsika::setup::Stack::ParticleType;
+using SetupTrack = corsika::setup::Trajectory;
namespace corsika::process::energy_loss {
@@ -53,7 +52,7 @@ namespace corsika::process::energy_loss {
* For shell correction, see Sec 6 of https://www.nap.edu/read/20066/chapter/8#115
*
*/
- HEPEnergyType EnergyLoss::BetheBloch(Particle& p, GrammageType const dX) {
+ HEPEnergyType EnergyLoss::BetheBloch(SetupParticle const& p, GrammageType const dX) {
// all these are material constants and have to come through Environment
// right now: values for nitrogen_D
@@ -70,20 +69,20 @@ namespace corsika::process::energy_loss {
// end of material constants
// this is the Bethe-Bloch coefficiet 4pi N_A r_e^2 m_e c^2
- auto const K = 0.307075_MeV / 1_mol * square(1_cm);
+ auto constexpr K = 0.307075_MeV / 1_mol * square(1_cm);
HEPEnergyType const E = p.GetEnergy();
HEPMassType const m = p.GetMass();
double const gamma = E / m;
- double const Z = p.GetChargeNumber();
- double const Z2 = pow(Z, 2);
- HEPMassType const me = particles::Electron::GetMass();
+ int const Z = p.GetChargeNumber();
+ int const Z2 = Z * Z;
+ HEPMassType constexpr me = particles::Electron::GetMass();
auto const m2 = m * m;
- auto const me2 = me * me;
- double const gamma2 = pow(gamma, 2);
+ auto constexpr me2 = me * me;
+ double const gamma2 = gamma * gamma;
double const beta2 = (gamma2 - 1) / gamma2; // 1-1/gamma2 (1-1/gamma)*(1+1/gamma);
// (gamma_2-1)/gamma_2 = (1-1/gamma2);
- double const c2 = 1; // HEP convention here c=c2=1
+ double constexpr c2 = 1; // HEP convention here c=c2=1
cout << "BetheBloch beta2=" << beta2 << " gamma2=" << gamma2 << endl;
[[maybe_unused]] double const eta2 = beta2 / (1 - beta2);
HEPMassType const Wmax =
@@ -148,7 +147,8 @@ namespace corsika::process::energy_loss {
(0.5 * log(aux) - beta2 - Cadj / Z - delta / 2 + barkas + bloch) * dX;
}
- process::EProcessReturn EnergyLoss::DoContinuous(Particle& p, Track& t) {
+ process::EProcessReturn EnergyLoss::DoContinuous(SetupParticle& p,
+ SetupTrack const& t) {
if (p.GetChargeNumber() == 0) return process::EProcessReturn::eOk;
GrammageType const dX =
p.GetNode()->GetModelProperties().IntegratedGrammage(t, t.GetLength());
@@ -170,12 +170,25 @@ namespace corsika::process::energy_loss {
p.SetEnergy(Enew);
MomentumUpdate(p, Enew);
fEnergyLossTot += dE;
- GetXbin(p, dE);
+ GetXbin(p, t, dE);
return status;
}
- units::si::LengthType EnergyLoss::MaxStepLength(Particle&, Track&) {
- return units::si::meter * std::numeric_limits<double>::infinity();
+ units::si::LengthType EnergyLoss::MaxStepLength(SetupParticle const& vParticle,
+ SetupTrack const& vTrack) const {
+ if (vParticle.GetChargeNumber() == 0) {
+ return units::si::meter * std::numeric_limits<double>::infinity();
+ }
+
+ auto constexpr dX = 1_g / square(1_cm);
+ auto const dE = -BetheBloch(vParticle, dX); // dE > 0
+ //~ auto const Ekin = vParticle.GetEnergy() - vParticle.GetMass();
+ auto const maxLoss = 0.01 * vParticle.GetEnergy();
+ auto const maxGrammage = maxLoss / dE * dX;
+
+ return vParticle.GetNode()->GetModelProperties().ArclengthFromGrammage(vTrack,
+ maxGrammage) *
+ 1.0001; // to make sure particle gets absorbed when DoContinuous() is called
}
void EnergyLoss::MomentumUpdate(corsika::setup::Stack::ParticleType& vP,
@@ -187,17 +200,17 @@ namespace corsika::process::energy_loss {
#include <corsika/geometry/CoordinateSystem.h>
- int EnergyLoss::GetXbin(corsika::setup::Stack::ParticleType& vP,
+ int EnergyLoss::GetXbin(SetupParticle const& vP, SetupTrack const& vTrack,
const HEPEnergyType dE) {
using namespace corsika::geometry;
CoordinateSystem const& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
- Point pos1(rootCS, 0_m, 0_m, 0_m);
- Point pos2(rootCS, 0_m, 0_m, vP.GetPosition().GetCoordinates()[2]);
- Vector delta = (pos2 - pos1) / 1_s;
- Trajectory t(Line(pos1, delta), 1_s);
+ Point const pos1(rootCS, 0_m, 0_m, 0_m);
+ Point const pos2(rootCS, 0_m, 0_m, vTrack.GetPosition(0).GetCoordinates()[2]);
+ auto const delta = (pos2 - pos1) / 1_s;
+ Trajectory const t(Line(pos1, delta), 1_s);
GrammageType const grammage =
vP.GetNode()->GetModelProperties().IntegratedGrammage(t, t.GetLength());
diff --git a/Processes/EnergyLoss/EnergyLoss.h b/Processes/EnergyLoss/EnergyLoss.h
index deaff4faf0837f6cac44ce2da975637393035d49..d0726d9d44220f27c356f5128a43a0092cecc558 100644
--- a/Processes/EnergyLoss/EnergyLoss.h
+++ b/Processes/EnergyLoss/EnergyLoss.h
@@ -25,34 +25,31 @@ namespace corsika::process::energy_loss {
class EnergyLoss : public corsika::process::ContinuousProcess<EnergyLoss> {
using MeVgcm2 = decltype(1e6 * units::si::electronvolt / units::si::gram *
- corsika::units::si::square(1e-2 * units::si::meter));
+ units::si::square(1e-2 * units::si::meter));
- void MomentumUpdate(corsika::setup::Stack::ParticleType&,
- corsika::units::si::HEPEnergyType Enew);
+ void MomentumUpdate(setup::Stack::ParticleType&, units::si::HEPEnergyType Enew);
public:
EnergyLoss();
void Init() {}
+ process::EProcessReturn DoContinuous(setup::Stack::ParticleType&,
+ setup::Trajectory const&);
+ units::si::LengthType MaxStepLength(setup::Stack::ParticleType const&,
+ setup::Trajectory const&) const;
- corsika::process::EProcessReturn DoContinuous(corsika::setup::Stack::ParticleType&,
- corsika::setup::Trajectory&);
- corsika::units::si::LengthType MaxStepLength(corsika::setup::Stack::ParticleType&,
- corsika::setup::Trajectory&);
-
- corsika::units::si::HEPEnergyType GetTotal() const { return fEnergyLossTot; }
+ units::si::HEPEnergyType GetTotal() const { return fEnergyLossTot; }
void PrintProfile() const;
private:
- corsika::units::si::HEPEnergyType BetheBloch(
- corsika::setup::Stack::ParticleType& p,
- const corsika::units::si::GrammageType dX);
+ static units::si::HEPEnergyType BetheBloch(setup::Stack::ParticleType const& p,
+ const units::si::GrammageType dX);
- int GetXbin(corsika::setup::Stack::ParticleType& p,
- const corsika::units::si::HEPEnergyType dE);
+ int GetXbin(setup::Stack::ParticleType const&, setup::Trajectory const&,
+ units::si::HEPEnergyType);
- corsika::units::si::HEPEnergyType fEnergyLossTot;
- corsika::units::si::GrammageType fdX; // profile binning
- std::map<int, double> fProfile; // longitudinal profile
+ units::si::HEPEnergyType fEnergyLossTot;
+ units::si::GrammageType fdX; // profile binning
+ std::map<int, double> fProfile; // longitudinal profile
};
} // namespace corsika::process::energy_loss
diff --git a/Processes/HadronicElasticModel/HadronicElasticModel.cc b/Processes/HadronicElasticModel/HadronicElasticModel.cc
index 33f162e41531bf4008f91b21bb98e8afc13a7abb..9e05c521320787b1f8a42ca3c154cfa708285769 100644
--- a/Processes/HadronicElasticModel/HadronicElasticModel.cc
+++ b/Processes/HadronicElasticModel/HadronicElasticModel.cc
@@ -18,14 +18,12 @@
#include <corsika/utl/COMBoost.h>
#include <corsika/setup/SetupStack.h>
-#include <corsika/setup/SetupTrajectory.h>
#include <iomanip>
#include <iostream>
using namespace corsika::setup;
-using Particle = corsika::setup::Stack::ParticleType;
-using Track = corsika::setup::Trajectory;
+using SetupParticle = corsika::setup::Stack::ParticleType;
namespace corsika::process::HadronicElasticModel {
@@ -38,7 +36,7 @@ namespace corsika::process::HadronicElasticModel {
template <>
units::si::GrammageType HadronicElasticInteraction::GetInteractionLength(
- Particle const& p, Track&) {
+ SetupParticle const& p) {
using namespace units::si;
if (p.GetPID() == particles::Code::Proton) {
auto const* currentNode = p.GetNode();
@@ -79,7 +77,7 @@ namespace corsika::process::HadronicElasticModel {
}
template <>
- process::EProcessReturn HadronicElasticInteraction::DoInteraction(Particle& p) {
+ process::EProcessReturn HadronicElasticInteraction::DoInteraction(SetupParticle& p) {
if (p.GetPID() != particles::Code::Proton) { return process::EProcessReturn::eOk; }
using namespace units::si;
diff --git a/Processes/HadronicElasticModel/HadronicElasticModel.h b/Processes/HadronicElasticModel/HadronicElasticModel.h
index 4180c04c514d00089fd044ab39027b9128a853ff..b3fd15fd2cdc6feac967182cff342b60f30e954d 100644
--- a/Processes/HadronicElasticModel/HadronicElasticModel.h
+++ b/Processes/HadronicElasticModel/HadronicElasticModel.h
@@ -56,8 +56,8 @@ namespace corsika::process::HadronicElasticModel {
units::si::CrossSectionType y = 0.05608 * units::si::barn);
void Init();
- template <typename Particle, typename Track>
- corsika::units::si::GrammageType GetInteractionLength(Particle const& p, Track&);
+ template <typename Particle>
+ corsika::units::si::GrammageType GetInteractionLength(Particle const& p);
template <typename Particle>
corsika::process::EProcessReturn DoInteraction(Particle&);
diff --git a/Processes/Pythia/Decay.cc b/Processes/Pythia/Decay.cc
index 222c68c8dd57d7fd57c453b62f53777e85ad062a..20d2aaa5044e484b2727e5cdc03f011fe36be83b 100644
--- a/Processes/Pythia/Decay.cc
+++ b/Processes/Pythia/Decay.cc
@@ -147,10 +147,6 @@ namespace corsika::process::pythia {
// set particle stable
Decay::SetStable(vP.GetPID());
-
- // remove original particle from corsika stack
- vP.Delete();
- // if (fCount>10) throw std::runtime_error("stop here");
}
} // namespace corsika::process::pythia
diff --git a/Processes/Pythia/Interaction.cc b/Processes/Pythia/Interaction.cc
index ae8104350086f5cfda94b174305d0f0366faa2ce..2a0cce386edc6dfc5a14f0f6fb84fd94aa5a85a7 100644
--- a/Processes/Pythia/Interaction.cc
+++ b/Processes/Pythia/Interaction.cc
@@ -15,7 +15,6 @@
#include <corsika/environment/NuclearComposition.h>
#include <corsika/geometry/FourVector.h>
#include <corsika/setup/SetupStack.h>
-#include <corsika/setup/SetupTrajectory.h>
#include <corsika/utl/COMBoost.h>
#include <tuple>
@@ -28,7 +27,6 @@ using namespace corsika;
using namespace corsika::setup;
using Projectile = corsika::setup::StackView::ParticleType;
using Particle = corsika::setup::Stack::ParticleType;
-using Track = Trajectory;
namespace corsika::process::pythia {
@@ -168,7 +166,7 @@ namespace corsika::process::pythia {
}
template <>
- units::si::GrammageType Interaction::GetInteractionLength(Particle& p, Track&) {
+ units::si::GrammageType Interaction::GetInteractionLength(Particle& p) {
using namespace units;
using namespace units::si;
@@ -411,8 +409,6 @@ namespace corsika::process::pythia {
<< "Elab_final=" << Elab_final / 1_GeV
<< ", Plab_final=" << (Plab_final / 1_GeV).GetComponents() << endl;
}
- // delete current particle
- vP.Delete();
}
return process::EProcessReturn::eOk;
}
diff --git a/Processes/Pythia/Interaction.h b/Processes/Pythia/Interaction.h
index ac27825acf0561d0d9e7e8e8b66961bdfc674e2f..3594358e0b39f528991183b4c414e8edaa3b5828 100644
--- a/Processes/Pythia/Interaction.h
+++ b/Processes/Pythia/Interaction.h
@@ -51,16 +51,16 @@ namespace corsika::process::pythia {
const corsika::particles::Code TargetId,
const corsika::units::si::HEPEnergyType CoMenergy);
- template <typename TParticle, typename TTrack>
- corsika::units::si::GrammageType GetInteractionLength(TParticle&, TTrack&);
+ template <typename TParticle>
+ corsika::units::si::GrammageType GetInteractionLength(TParticle&);
/**
In this function PYTHIA is called to produce one event. The
event is copied (and boosted) into the shower lab frame.
*/
- template <typename TProjctile>
- corsika::process::EProcessReturn DoInteraction(TProjctile&);
+ template <typename TProjectile>
+ corsika::process::EProcessReturn DoInteraction(TProjectile&);
private:
corsika::random::RNG& fRNG =
diff --git a/Processes/Pythia/Random.cc b/Processes/Pythia/Random.cc
index 8fcbcf4ee7274ba36252f7969306073f3cffddd5..0d08403ca1d890a062e2632124079c6a3cb220f6 100644
--- a/Processes/Pythia/Random.cc
+++ b/Processes/Pythia/Random.cc
@@ -12,9 +12,6 @@
namespace corsika::process::pythia {
- double Random::flat() {
- std::uniform_real_distribution<double> dist;
- return dist(fRNG);
- }
+ double Random::flat() { return fDist(fRNG); }
} // namespace corsika::process::pythia
diff --git a/Processes/Pythia/Random.h b/Processes/Pythia/Random.h
index 276fc532ca53ac84b708e8ad12cb16fbaed01c9a..cd35fc4eeafcbdf4ec507cf54ab966271a8183f0 100644
--- a/Processes/Pythia/Random.h
+++ b/Processes/Pythia/Random.h
@@ -23,6 +23,7 @@ namespace corsika::process {
double flat();
private:
+ std::uniform_real_distribution<double> fDist;
corsika::random::RNG& fRNG =
corsika::random::RNGManager::GetInstance().GetRandomStream("pythia");
};
diff --git a/Processes/Pythia/testPythia.cc b/Processes/Pythia/testPythia.cc
index 3fb2df13ecd1c24e674997a3c3aa468164f21d06..f6f09b84c17f80f8b579d1a11eba843e6ddbc964 100644
--- a/Processes/Pythia/testPythia.cc
+++ b/Processes/Pythia/testPythia.cc
@@ -118,14 +118,6 @@ TEST_CASE("pythia process") {
auto const* nodePtr = theMedium.get(); // save the medium for later use before moving it
universe.AddChild(std::move(theMedium));
- geometry::Point const origin(cs, {0_m, 0_m, 0_m});
- geometry::Vector<units::si::SpeedType::dimension_type> v(cs, 0_m / second, 0_m / second,
- 1_m / second);
- geometry::Line line(origin, v);
- geometry::Trajectory<geometry::Line> track(line, 10_s);
-
- random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
-
SECTION("pythia decay") {
setup::Stack stack;
@@ -173,7 +165,6 @@ TEST_CASE("pythia process") {
process::pythia::Interaction model;
model.Init();
[[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(projectile);
- [[maybe_unused]] const GrammageType length =
- model.GetInteractionLength(particle, track);
+ [[maybe_unused]] const GrammageType length = model.GetInteractionLength(particle);
}
}
diff --git a/Processes/Sibyll/Interaction.cc b/Processes/Sibyll/Interaction.cc
index d674f9a3a775970bf60de08c0cf33b55fc73a17d..c089ec083c9c7d9ec1d91136cb0c7796dc958c57 100644
--- a/Processes/Sibyll/Interaction.cc
+++ b/Processes/Sibyll/Interaction.cc
@@ -118,8 +118,8 @@ namespace corsika::process::sibyll {
}
template <>
- units::si::GrammageType Interaction::GetInteractionLength(SetupParticle& vP,
- Track&) const {
+ units::si::GrammageType Interaction::GetInteractionLength(
+ SetupParticle const& vP) const {
using namespace units;
using namespace units::si;
@@ -366,8 +366,6 @@ namespace corsika::process::sibyll {
<< ", Plab_final=" << (Plab_final / 1_GeV).GetComponents() << endl;
}
}
- // delete current particle
- vP.Delete();
return process::EProcessReturn::eOk;
}
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
index e7d7a0902116c2fceb9240cf9266d04a2ea2fd72..643bd0d51bd4f387aafa3870f2dddfb0ed1fe9d5 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/Processes/Sibyll/Interaction.h
@@ -52,8 +52,8 @@ namespace corsika::process::sibyll {
GetCrossSection(const corsika::particles::Code, const corsika::particles::Code,
const corsika::units::si::HEPEnergyType) const;
- template <typename TParticle, typename TTrack>
- corsika::units::si::GrammageType GetInteractionLength(TParticle&, TTrack&) const;
+ template <typename TParticle>
+ corsika::units::si::GrammageType GetInteractionLength(TParticle const&) const;
/**
In this function SIBYLL is called to produce one event. The
diff --git a/Processes/Sibyll/NuclearInteraction.cc b/Processes/Sibyll/NuclearInteraction.cc
index 31abd9da350e1c80c24868eb2313e860b51bd821..907e95b271bb44cfb4d2982f0f1f40f88cbb0eca 100644
--- a/Processes/Sibyll/NuclearInteraction.cc
+++ b/Processes/Sibyll/NuclearInteraction.cc
@@ -217,7 +217,7 @@ namespace corsika::process::sibyll {
template <>
template <>
units::si::GrammageType NuclearInteraction<SetupEnvironment>::GetInteractionLength(
- Particle& vP, Track&) {
+ Particle& vP) {
using namespace units;
using namespace units::si;
diff --git a/Processes/Sibyll/NuclearInteraction.h b/Processes/Sibyll/NuclearInteraction.h
index 6bee9ba7a7964e74112c2bf60e5cbd4b1c23c2fa..261f931321ed98119ab02f4e3cace952308fdbb7 100644
--- a/Processes/Sibyll/NuclearInteraction.h
+++ b/Processes/Sibyll/NuclearInteraction.h
@@ -53,11 +53,11 @@ namespace corsika::process::sibyll {
std::tuple<corsika::units::si::CrossSectionType, corsika::units::si::CrossSectionType>
GetCrossSection(Particle& p, const corsika::particles::Code TargetId);
- template <typename Particle, typename Track>
- corsika::units::si::GrammageType GetInteractionLength(Particle&, Track&);
+ template <typename Particle>
+ corsika::units::si::GrammageType GetInteractionLength(Particle&);
- template <typename Projectle>
- corsika::process::EProcessReturn DoInteraction(Projectle&);
+ template <typename Projectile>
+ corsika::process::EProcessReturn DoInteraction(Projectile&);
private:
TEnvironment const& fEnvironment;
diff --git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc
index d5f5d9bc8e39ac775de3fb9d541a18b596941bad..0c42c90abb86a758c45f5ad4a79ebe87d3f1d3fa 100644
--- a/Processes/Sibyll/testSibyll.cc
+++ b/Processes/Sibyll/testSibyll.cc
@@ -98,12 +98,6 @@ TEST_CASE("SibyllInterface", "[processes]") {
const geometry::CoordinateSystem& cs = env.GetCoordinateSystem();
- geometry::Point const origin(cs, {0_m, 0_m, 0_m});
- geometry::Vector<units::si::SpeedType::dimension_type> v(cs, 0_m / second, 0_m / second,
- 1_m / second);
- geometry::Line line(origin, v);
- geometry::Trajectory<geometry::Line> track(line, 10_s);
-
random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
SECTION("InteractionInterface") {
@@ -126,8 +120,7 @@ TEST_CASE("SibyllInterface", "[processes]") {
model.Init();
[[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(projectile);
- [[maybe_unused]] const GrammageType length =
- model.GetInteractionLength(particle, track);
+ [[maybe_unused]] const GrammageType length = model.GetInteractionLength(particle);
}
SECTION("NuclearInteractionInterface") {
@@ -153,8 +146,7 @@ TEST_CASE("SibyllInterface", "[processes]") {
model.Init();
[[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(projectile);
- [[maybe_unused]] const GrammageType length =
- model.GetInteractionLength(particle, track);
+ [[maybe_unused]] const GrammageType length = model.GetInteractionLength(particle);
}
SECTION("DecayInterface") {
diff --git a/Processes/StackInspector/StackInspector.cc b/Processes/StackInspector/StackInspector.cc
index 814763d5cbfcbc0bd4201db70f220dc4e8ea123f..e56503378ed5abfc15578c474ed0afd4d54d6412 100644
--- a/Processes/StackInspector/StackInspector.cc
+++ b/Processes/StackInspector/StackInspector.cc
@@ -37,7 +37,7 @@ template <typename TStack>
StackInspector<TStack>::~StackInspector() {}
template <typename TStack>
-process::EProcessReturn StackInspector<TStack>::DoStack(TStack& vS) {
+process::EProcessReturn StackInspector<TStack>::DoStack(TStack const& vS) {
if (!fReport) return process::EProcessReturn::eOk;
[[maybe_unused]] int i = 0;
HEPEnergyType Etot = 0_GeV;
diff --git a/Processes/StackInspector/StackInspector.h b/Processes/StackInspector/StackInspector.h
index 57a929dd2362b302d43d05175878a257c5c8a832..62f36eeebae774f74d26953697707e698c75da6e 100644
--- a/Processes/StackInspector/StackInspector.h
+++ b/Processes/StackInspector/StackInspector.h
@@ -30,7 +30,7 @@ namespace corsika::process {
~StackInspector();
void Init();
- EProcessReturn DoStack(TStack&);
+ EProcessReturn DoStack(TStack const&);
private:
bool fReport;
diff --git a/Processes/SwitchProcess/CMakeLists.txt b/Processes/SwitchProcess/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..feac9f360ced5356d90ecd87cc7d4ca51939e0f6
--- /dev/null
+++ b/Processes/SwitchProcess/CMakeLists.txt
@@ -0,0 +1,50 @@
+set (
+ MODEL_HEADERS
+ SwitchProcess.h
+ )
+
+set (
+ MODEL_NAMESPACE
+ corsika/process/switch_process
+ )
+
+add_library (ProcessSwitch INTERFACE)
+CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessSwitch ${MODEL_NAMESPACE} ${MODEL_HEADERS})
+
+set_target_properties (
+ ProcessStackInspector
+ PROPERTIES
+ VERSION ${PROJECT_VERSION}
+ SOVERSION 1
+# PUBLIC_HEADER "${MODEL_HEADERS}"
+ )
+
+# target dependencies on other libraries (also the header onlys)
+target_link_libraries (
+ ProcessSwitch
+ INTERFACE
+ CORSIKAunits
+ CORSIKAprocesssequence
+ )
+
+target_include_directories (
+ ProcessSwitch
+ INTERFACE
+ $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
+ $<INSTALL_INTERFACE:include/include>
+ )
+
+install (FILES ${MODEL_HEADERS} DESTINATION include/${MODEL_NAMESPACE})
+
+# --------------------
+# code unit testing
+add_executable (testSwitchProcess testSwitchProcess.cc)
+
+target_link_libraries (
+ testSwitchProcess
+ ProcessSwitch
+ CORSIKAstackinterface
+ CORSIKAthirdparty # for catch2
+)
+
+CORSIKA_ADD_TEST(testSwitchProcess)
diff --git a/Processes/SwitchProcess/SwitchProcess.h b/Processes/SwitchProcess/SwitchProcess.h
new file mode 100644
index 0000000000000000000000000000000000000000..31e59599bc6dceb68edf4ad030d284eb37c69bce
--- /dev/null
+++ b/Processes/SwitchProcess/SwitchProcess.h
@@ -0,0 +1,98 @@
+#ifndef _corsika_SwitchProcess_h
+#define _corsika_SwitchProcess_h
+
+#include <corsika/process/InteractionProcess.h>
+#include <corsika/process/ProcessSequence.h>
+#include <corsika/setup/SetupStack.h>
+#include <corsika/units/PhysicalUnits.h>
+
+namespace corsika::process::switch_process {
+
+ /**
+ * This process provides an energy-based switch between two interaction processes P1 and
+ * P1. For energies below the threshold, P1 is invoked, otherwise P2. Both can be either
+ * single interaction processes or multiple ones combined in a ProcessSequence. A
+ * SwitchProcess itself will always be regarded as a ProcessSequence rather than an
+ * InteractionProcess when assembled into a greater ProcessSequence.
+ */
+
+ template <class TLowEProcess, class THighEProcess>
+ class SwitchProcess : public BaseProcess<SwitchProcess<TLowEProcess, THighEProcess>> {
+ TLowEProcess& fLowEProcess;
+ THighEProcess& fHighEProcess;
+ units::si::HEPEnergyType const fThresholdEnergy;
+
+ public:
+ SwitchProcess(TLowEProcess& vLowEProcess, THighEProcess& vHighEProcess,
+ units::si::HEPEnergyType vThresholdEnergy)
+ : fLowEProcess(vLowEProcess)
+ , fHighEProcess(vHighEProcess)
+ , fThresholdEnergy(vThresholdEnergy) {}
+
+ void Init() {
+ fLowEProcess.Init();
+ fHighEProcess.Init();
+ }
+
+ template <typename TParticle>
+ corsika::units::si::InverseGrammageType GetInverseInteractionLength(TParticle& p) {
+ return 1 / GetInteractionLength(p);
+ }
+
+ template <typename TParticle>
+ units::si::GrammageType GetInteractionLength(TParticle& vParticle) {
+ if (vParticle.GetEnergy() < fThresholdEnergy) {
+ if constexpr (is_process_sequence_v<TLowEProcess>) {
+ return fLowEProcess.GetTotalInteractionLength(vParticle);
+ } else {
+ return fLowEProcess.GetInteractionLength(vParticle);
+ }
+ } else {
+ if constexpr (is_process_sequence_v<THighEProcess>) {
+ return fHighEProcess.GetTotalInteractionLength(vParticle);
+ } else {
+ return fHighEProcess.GetInteractionLength(vParticle);
+ }
+ }
+ }
+
+ // required to conform to ProcessSequence interface. We cannot just
+ // implement DoInteraction() because we want to call SelectInteraction
+ // in case a member process is a ProcessSequence.
+ template <typename TParticle, typename TSecondaries>
+ EProcessReturn SelectInteraction(
+ TParticle& vP, TSecondaries& vS,
+ [[maybe_unused]] corsika::units::si::InverseGrammageType lambda_select,
+ corsika::units::si::InverseGrammageType& lambda_inv_count) {
+ if (vP.GetEnergy() < fThresholdEnergy) {
+ if constexpr (is_process_sequence_v<TLowEProcess>) {
+ return fLowEProcess.SelectInteraction(vP, vS, lambda_select, lambda_inv_count);
+ } else {
+ lambda_inv_count += fLowEProcess.GetInverseInteractionLength(vP);
+ // check if we should execute THIS process and then EXIT
+ if (lambda_select < lambda_inv_count) {
+ fLowEProcess.DoInteraction(vS);
+ return EProcessReturn::eInteracted;
+ } else {
+ return EProcessReturn::eOk;
+ }
+ }
+ } else {
+ if constexpr (is_process_sequence_v<THighEProcess>) {
+ return fHighEProcess.SelectInteraction(vP, vS, lambda_select, lambda_inv_count);
+ } else {
+ lambda_inv_count += fHighEProcess.GetInverseInteractionLength(vP);
+ // check if we should execute THIS process and then EXIT
+ if (lambda_select < lambda_inv_count) {
+ fHighEProcess.DoInteraction(vS);
+ return EProcessReturn::eInteracted;
+ } else {
+ return EProcessReturn::eOk;
+ }
+ }
+ }
+ }
+ };
+} // namespace corsika::process::switch_process
+
+#endif
diff --git a/Processes/SwitchProcess/testSwitchProcess.cc b/Processes/SwitchProcess/testSwitchProcess.cc
new file mode 100644
index 0000000000000000000000000000000000000000..4e902e8857367bfbb4d1ab7643941cdcce27c04c
--- /dev/null
+++ b/Processes/SwitchProcess/testSwitchProcess.cc
@@ -0,0 +1,256 @@
+#include <corsika/process/switch_process/SwitchProcess.h>
+#include <corsika/stack/SecondaryView.h>
+#include <corsika/stack/Stack.h>
+#include <corsika/units/PhysicalUnits.h>
+
+#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
+ // cpp file
+#include <catch2/catch.hpp>
+
+#include <algorithm>
+#include <random>
+
+using namespace corsika;
+using namespace corsika::process;
+using namespace corsika::units::si;
+
+class TestStackData {
+
+public:
+ // these functions are needed for the Stack interface
+ void Init() {}
+ void Clear() { fData.clear(); }
+ unsigned int GetSize() const { return fData.size(); }
+ unsigned int GetCapacity() const { return fData.size(); }
+ void Copy(int i1, int i2) { fData[i2] = fData[i1]; }
+ void Swap(int i1, int i2) { std::swap(fData[i1], fData[i2]); }
+
+ // custom data access function
+ void SetData(unsigned int i, HEPEnergyType v) { fData[i] = v; }
+ HEPEnergyType GetData(const unsigned int i) const { return fData[i]; }
+
+ // these functions are also needed by the Stack interface
+ void IncrementSize() { fData.resize(fData.size() + 1); }
+ void DecrementSize() {
+ if (fData.size() > 0) { fData.pop_back(); }
+ }
+
+ // custom private data section
+private:
+ std::vector<HEPEnergyType> fData;
+};
+
+/**
+ * From static_cast of a StackIterator over entries in the
+ * TestStackData class you get and object of type
+ * TestParticleInterface defined here
+ *
+ * It provides Get/Set methods to read and write data to the "Data"
+ * storage of TestStackData obtained via
+ * "StackIteratorInterface::GetStackData()", given the index of the
+ * iterator "StackIteratorInterface::GetIndex()"
+ *
+ */
+template <typename StackIteratorInterface>
+class TestParticleInterface
+ : public corsika::stack::ParticleBase<StackIteratorInterface> {
+
+public:
+ using corsika::stack::ParticleBase<StackIteratorInterface>::GetStackData;
+ using corsika::stack::ParticleBase<StackIteratorInterface>::GetIndex;
+
+ /*
+ The SetParticleData methods are called for creating new entries
+ on the stack. You can specifiy various parametric versions to
+ perform this task:
+ */
+
+ // default version for particle-creation from input data
+ void SetParticleData(const std::tuple<HEPEnergyType> v) { SetEnergy(std::get<0>(v)); }
+ void SetParticleData(TestParticleInterface<StackIteratorInterface>& /*parent*/,
+ std::tuple<HEPEnergyType> v) {
+ SetEnergy(std::get<0>(v));
+ }
+
+ // here are the fundamental methods for access to TestStackData data
+ void SetEnergy(HEPEnergyType v) { GetStackData().SetData(GetIndex(), v); }
+ HEPEnergyType GetEnergy() const { return GetStackData().GetData(GetIndex()); }
+};
+
+using SimpleStack = corsika::stack::Stack<TestStackData, TestParticleInterface>;
+
+// see issue 161
+#if defined(__clang__)
+using StackTestView = corsika::stack::SecondaryView<TestStackData, TestParticleInterface>;
+#elif defined(__GNUC__) || defined(__GNUG__)
+using StackTestView = corsika::stack::MakeView<SimpleStack>::type;
+#endif
+
+auto constexpr kgMSq = 1_kg / (1_m * 1_m);
+
+template <int N>
+struct DummyProcess : InteractionProcess<DummyProcess<N>> {
+ void Init() {}
+
+ template <typename TParticle>
+ corsika::units::si::GrammageType GetInteractionLength(TParticle const&) const {
+ return N * kgMSq;
+ }
+
+ template <typename TSecondaries>
+ corsika::process::EProcessReturn DoInteraction(TSecondaries& vSec) {
+ // to figure out which process was selected in the end, we produce N
+ // secondaries for DummyProcess<N>
+
+ for (int i = 0; i < N; ++i) {
+ vSec.AddSecondary(std::tuple<HEPEnergyType>{vSec.GetEnergy() / N});
+ }
+
+ return EProcessReturn::eOk;
+ }
+};
+
+using DummyLowEnergyProcess = DummyProcess<1>;
+using DummyHighEnergyProcess = DummyProcess<2>;
+using DummyAdditionalProcess = DummyProcess<3>;
+
+TEST_CASE("SwitchProcess from InteractionProcess") {
+ DummyLowEnergyProcess low;
+ DummyHighEnergyProcess high;
+ DummyAdditionalProcess proc;
+
+ switch_process::SwitchProcess switchProcess(low, high, 1_TeV);
+ auto seq = switchProcess << proc;
+
+ SimpleStack stack;
+
+ SECTION("low energy") {
+ stack.AddParticle(std::tuple<HEPEnergyType>{0.5_TeV});
+ auto p = stack.GetNextParticle();
+
+ // low energy process returns 1 kg/m²
+ SECTION("interaction length") {
+ REQUIRE(switchProcess.GetInteractionLength(p) / kgMSq == Approx(1));
+ REQUIRE(seq.GetTotalInteractionLength(p) / kgMSq == Approx(3. / 4));
+ }
+
+ // low energy process creates 1 secondary
+ //~ SECTION("SelectInteraction") {
+ //~ typename SimpleStack::ParticleType theParticle =
+ //~ stack.GetNextParticle(); // as in corsika::Cascade
+ //~ StackTestView view(theParticle);
+ //~ auto projectile = view.GetProjectile();
+
+ //~ InverseGrammageType invLambda = 0 / kgMSq;
+ //~ switchProcess.SelectInteraction(p, projectile, 0.01 / kgMSq, invLambda);
+
+ //~ REQUIRE(view.GetSize() == 1);
+ //~ }
+ }
+
+ SECTION("high energy") {
+ stack.AddParticle(std::tuple<HEPEnergyType>{4_TeV});
+ auto p = stack.GetNextParticle();
+
+ // high energy process returns 2 kg/m²
+ SECTION("interaction length") {
+ REQUIRE(switchProcess.GetInteractionLength(p) / kgMSq == Approx(2));
+ REQUIRE(seq.GetTotalInteractionLength(p) / kgMSq == Approx(6. / 5));
+ }
+
+ // high energy process creates 2 secondaries
+ SECTION("SelectInteraction") {
+ typename SimpleStack::ParticleType theParticle =
+ stack.GetNextParticle(); // as in corsika::Cascade
+ StackTestView view(theParticle);
+ auto projectile = view.GetProjectile();
+
+ InverseGrammageType invLambda = 0 / kgMSq;
+ switchProcess.SelectInteraction(p, projectile, 0.01 / kgMSq, invLambda);
+
+ REQUIRE(view.GetSize() == 2);
+ }
+ }
+}
+
+TEST_CASE("SwitchProcess from ProcessSequence") {
+ DummyProcess<1> innerA;
+ DummyProcess<2> innerB;
+ DummyProcess<3> outer;
+ DummyProcess<4> additional;
+
+ auto seq = innerA << innerB;
+
+ switch_process::SwitchProcess switchProcess(seq, outer, 1_TeV);
+ auto completeSeq = switchProcess << additional;
+
+ SimpleStack stack;
+
+ SECTION("low energy") {
+ stack.AddParticle(std::tuple<HEPEnergyType>{0.5_TeV});
+ auto p = stack.GetNextParticle();
+
+ SECTION("interaction length") {
+ REQUIRE(switchProcess.GetInteractionLength(p) / kgMSq == Approx(2. / 3));
+ REQUIRE(completeSeq.GetTotalInteractionLength(p) / kgMSq == Approx(4. / 7));
+ }
+
+ SECTION("SelectInteraction") {
+ std::vector<int> numberOfSecondaries;
+
+ for (int i = 0; i < 1000; ++i) {
+ typename SimpleStack::ParticleType theParticle =
+ stack.GetNextParticle(); // as in corsika::Cascade
+ StackTestView view(theParticle);
+ auto projectile = view.GetProjectile();
+
+ double r = i / 1000.;
+ InverseGrammageType invLambda = r * 7./4 / kgMSq;
+
+ InverseGrammageType accumulator = 0 / kgMSq;
+ completeSeq.SelectInteraction(p, projectile, invLambda, accumulator);
+
+ numberOfSecondaries.push_back(view.GetSize());
+ }
+
+ auto const mean =
+ std::accumulate(numberOfSecondaries.cbegin(), numberOfSecondaries.cend(), 0.) /
+ numberOfSecondaries.size();
+ REQUIRE(mean == Approx(12./7.).margin(0.01));
+ }
+ }
+
+ SECTION("high energy") {
+ stack.AddParticle(std::tuple<HEPEnergyType>{3.0_TeV});
+ auto p = stack.GetNextParticle();
+
+ SECTION("interaction length") {
+ REQUIRE(switchProcess.GetInteractionLength(p) / kgMSq == Approx(3));
+ REQUIRE(completeSeq.GetTotalInteractionLength(p) / kgMSq == Approx(12. / 7.));
+ }
+
+ SECTION("SelectInteraction") {
+ std::vector<int> numberOfSecondaries;
+
+ for (int i = 0; i < 1000; ++i) {
+ typename SimpleStack::ParticleType theParticle =
+ stack.GetNextParticle(); // as in corsika::Cascade
+ StackTestView view(theParticle);
+ auto projectile = view.GetProjectile();
+
+ double r = i / 1000.;
+ InverseGrammageType invLambda = r * 7. / 12. / kgMSq;
+
+ InverseGrammageType accumulator = 0 / kgMSq;
+ completeSeq.SelectInteraction(p, projectile, invLambda, accumulator);
+
+ numberOfSecondaries.push_back(view.GetSize());
+ }
+
+ auto const mean =
+ std::accumulate(numberOfSecondaries.cbegin(), numberOfSecondaries.cend(), 0.) /
+ numberOfSecondaries.size();
+ REQUIRE(mean == Approx(24. / 7.).margin(0.01));
+ }
+ }
+}
diff --git a/Processes/TrackWriter/TrackWriter.cc b/Processes/TrackWriter/TrackWriter.cc
index e422bce8af8f054819cd03c174eaf594e85d6c76..d3247e87c8ad0e66620afcad2b2dc18d3e1c6049 100644
--- a/Processes/TrackWriter/TrackWriter.cc
+++ b/Processes/TrackWriter/TrackWriter.cc
@@ -37,9 +37,9 @@ namespace corsika::process::track_writer {
using namespace units::si;
auto const start = vT.GetPosition(0).GetCoordinates();
auto const delta = vT.GetPosition(1).GetCoordinates() - start;
- auto const& name = particles::GetName(vP.GetPID());
+ auto const pdg = static_cast<int>(particles::GetPDG(vP.GetPID()));
- fFile << name << " " << vP.GetEnergy() / 1_eV << ' ' << start[0] / 1_m << ' '
+ fFile << pdg << ' ' << vP.GetEnergy() / 1_eV << ' ' << start[0] / 1_m << ' '
<< start[1] / 1_m << ' ' << start[2] / 1_m << " " << delta[0] / 1_m << ' '
<< delta[1] / 1_m << ' ' << delta[2] / 1_m << '\n';
diff --git a/Processes/TrackingLine/TrackingLine.cc b/Processes/TrackingLine/TrackingLine.cc
index 319030f91d5347d47af58cddbdea607684a16677..149556bf245b9fb418fb3db03f267e64b688fb5f 100644
--- a/Processes/TrackingLine/TrackingLine.cc
+++ b/Processes/TrackingLine/TrackingLine.cc
@@ -16,31 +16,27 @@
#include <corsika/geometry/Vector.h>
#include <corsika/process/tracking_line/TrackingLine.h>
-#include <algorithm>
-#include <iostream>
+#include <limits>
#include <stdexcept>
#include <utility>
-using namespace corsika;
+using namespace corsika::geometry;
+using namespace corsika::units::si;
namespace corsika::process::tracking_line {
- std::optional<std::pair<corsika::units::si::TimeType, corsika::units::si::TimeType>>
- TimeOfIntersection(corsika::geometry::Line const& line,
- geometry::Sphere const& sphere) {
+ std::optional<std::pair<TimeType, TimeType>> TimeOfIntersection(Line const& line,
+ Sphere const& sphere) {
auto const delta = line.GetR0() - sphere.GetCenter();
auto const v = line.GetV0();
- auto const vSqNorm = v.squaredNorm();
+ auto const vSqNorm =
+ v.squaredNorm(); // todo: get rid of this by having V0 normalized always
auto const R = sphere.GetRadius();
auto const vDotDelta = v.dot(delta);
auto const discriminant =
vDotDelta * vDotDelta - vSqNorm * (delta.squaredNorm() - R * R);
- //~ std::cout << "discriminant: " << discriminant << std::endl;
- //~ std::cout << "alpha: " << alpha << std::endl;
- //~ std::cout << "beta: " << beta << std::endl;
-
if (discriminant.magnitude() > 0) {
auto const sqDisc = sqrt(discriminant);
auto const invDenom = 1 / vSqNorm;
@@ -50,4 +46,17 @@ namespace corsika::process::tracking_line {
return {};
}
}
+
+ TimeType TimeOfIntersection(Line const& vLine, Plane const& vPlane) {
+ auto const delta = vPlane.GetCenter() - vLine.GetR0();
+ auto const v = vLine.GetV0();
+ auto const n = vPlane.GetNormal();
+ auto const c = n.dot(v);
+
+ if (c.magnitude() == 0) {
+ return std::numeric_limits<TimeType::value_type>::infinity() * 1_s;
+ } else {
+ return n.dot(delta) / c;
+ }
+ }
} // namespace corsika::process::tracking_line
diff --git a/Processes/TrackingLine/TrackingLine.h b/Processes/TrackingLine/TrackingLine.h
index a801a2ee2baa206424211151c6a7fb9a616d0b70..df02264fe6f4c5640fe40c799255e87ce26a958e 100644
--- a/Processes/TrackingLine/TrackingLine.h
+++ b/Processes/TrackingLine/TrackingLine.h
@@ -13,6 +13,7 @@
#define _include_corsika_processes_TrackingLine_h_
#include <corsika/geometry/Line.h>
+#include <corsika/geometry/Plane.h>
#include <corsika/geometry/Sphere.h>
#include <corsika/geometry/Trajectory.h>
#include <corsika/geometry/Vector.h>
@@ -33,8 +34,10 @@ namespace corsika::process {
namespace tracking_line {
std::optional<std::pair<corsika::units::si::TimeType, corsika::units::si::TimeType>>
- TimeOfIntersection(corsika::geometry::Line const& line,
- geometry::Sphere const& sphere);
+ TimeOfIntersection(geometry::Line const&, geometry::Sphere const&);
+
+ corsika::units::si::TimeType TimeOfIntersection(geometry::Line const&,
+ geometry::Plane const&);
class TrackingLine {
diff --git a/Processes/UrQMD/UrQMD.cc b/Processes/UrQMD/UrQMD.cc
index df4456bfb965f398720e33aafd286ed56bbdb307..98f7a7385727d30836e7590bfb6bd54f8de4f542 100644
--- a/Processes/UrQMD/UrQMD.cc
+++ b/Processes/UrQMD/UrQMD.cc
@@ -2,9 +2,6 @@
#include <corsika/geometry/Vector.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/process/urqmd/UrQMD.h>
-#include <corsika/random/RNGManager.h>
-#include <corsika/setup/SetupStack.h>
-#include <corsika/setup/SetupTrajectory.h>
#include <corsika/units/PhysicalUnits.h>
#include <algorithm>
@@ -20,7 +17,6 @@ UrQMD::UrQMD() { iniurqmd_(); }
using SetupStack = corsika::setup::Stack;
using SetupParticle = corsika::setup::Stack::StackIterator;
using SetupProjectile = corsika::setup::StackView::StackIterator;
-using SetupTrack = corsika::setup::Trajectory;
template <typename TParticle> // need template here, as this is called both with
// SetupParticle as well as SetupProjectile
@@ -82,7 +78,7 @@ bool UrQMD::CanInteract(particles::Code vCode) const {
vCode) != std::cend(validProjectileCodes);
}
-GrammageType UrQMD::GetInteractionLength(SetupParticle& vParticle, SetupTrack&) const {
+GrammageType UrQMD::GetInteractionLength(SetupParticle& vParticle) const {
if (!CanInteract(vParticle.GetPID())) {
// we could do the canInteract check in GetCrossSection, too but if
// we do it here we have the advantage of avoiding the loop
diff --git a/Processes/UrQMD/UrQMD.h b/Processes/UrQMD/UrQMD.h
index 0b34e2e147a94c58a933d0cf508a428be5c7b917..a77e6e843b38055bc96c470c79f2902fcfa58ce5 100644
--- a/Processes/UrQMD/UrQMD.h
+++ b/Processes/UrQMD/UrQMD.h
@@ -5,7 +5,6 @@
#include <corsika/process/InteractionProcess.h>
#include <corsika/random/RNGManager.h>
#include <corsika/setup/SetupStack.h>
-#include <corsika/setup/SetupTrajectory.h>
#include <corsika/units/PhysicalUnits.h>
#include <array>
@@ -15,8 +14,9 @@ namespace corsika::process::UrQMD {
class UrQMD : public corsika::process::InteractionProcess<UrQMD> {
public:
UrQMD();
+ void Init() {}
corsika::units::si::GrammageType GetInteractionLength(
- corsika::setup::Stack::StackIterator&, corsika::setup::Trajectory&) const;
+ corsika::setup::Stack::StackIterator&) const;
template <typename TParticle>
corsika::units::si::CrossSectionType GetCrossSection(TParticle const&,
diff --git a/README.md b/README.md
index 67dd505f0d4dafc3f5ee73e055df7e4754429c99..d5cf8f7c2764daf2fe25f5c861502318a7ea4026 100644
--- a/README.md
+++ b/README.md
@@ -60,7 +60,7 @@ CORSIKA 8 is tested regularly at least on gcc7.3.0 and clang-6.0.0.
Additional software prerequisites: eigen3, boost, cmake, g++, git.
On a bare Ubuntu 18.04, just add:
```
-sudo apt-get install libeigen3-dev libboost-dev cmake g++ git
+sudo apt install libeigen3-dev cmake g++ git
```
Follow these steps to download and install CORSIKA 8 milestone2