diff --git a/Environment/Environment.h b/Environment/Environment.h
index 335f776ee69d76bc77fb2fef628761ae8b67508c..1aea20cda07f642cabb82a0fb9eb18b6a73d0e4a 100644
--- a/Environment/Environment.h
+++ b/Environment/Environment.h
@@ -14,8 +14,8 @@
#include <corsika/environment/IMediumModel.h>
#include <corsika/environment/VolumeTreeNode.h>
#include <corsika/geometry/Point.h>
-#include <corsika/geometry/Sphere.h>
#include <corsika/geometry/RootCoordinateSystem.h>
+#include <corsika/geometry/Sphere.h>
#include <corsika/setup/SetupEnvironment.h>
#include <limits>
@@ -23,9 +23,9 @@ namespace corsika::environment {
struct Universe : public corsika::geometry::Sphere {
Universe(corsika::geometry::CoordinateSystem const& pCS)
: corsika::geometry::Sphere(
- corsika::geometry::Point{
- pCS, 0 * corsika::units::si::meter,
- 0 * corsika::units::si::meter, 0 * corsika::units::si::meter},
+ corsika::geometry::Point{pCS, 0 * corsika::units::si::meter,
+ 0 * corsika::units::si::meter,
+ 0 * corsika::units::si::meter},
corsika::units::si::meter * std::numeric_limits<double>::infinity()) {}
bool Contains(corsika::geometry::Point const&) const override { return true; }
@@ -35,8 +35,9 @@ namespace corsika::environment {
class Environment {
public:
Environment()
- : fCoordinateSystem{corsika::geometry::RootCoordinateSystem::GetInstance().GetRootCS()},
- fUniverse(std::make_unique<VolumeTreeNode<IEnvironmentModel>>(
+ : fCoordinateSystem{corsika::geometry::RootCoordinateSystem::GetInstance()
+ .GetRootCS()}
+ , fUniverse(std::make_unique<VolumeTreeNode<IEnvironmentModel>>(
std::make_unique<Universe>(fCoordinateSystem))) {}
using IEnvironmentModel = corsika::setup::IEnvironmentModel;
diff --git a/Environment/HomogeneousMedium.h b/Environment/HomogeneousMedium.h
index 019774d5120a605fade3c50c62c7baa88bb8eb00..566b565e7ed4b7d03f37cf9a4a7c0baf7b232409 100644
--- a/Environment/HomogeneousMedium.h
+++ b/Environment/HomogeneousMedium.h
@@ -45,10 +45,10 @@ namespace corsika::environment {
corsika::geometry::Trajectory<corsika::geometry::Line> const& pTraj,
corsika::units::si::LengthType pTo) const override {
using namespace corsika::units::si;
- pTo * fDensity;
+ return pTo * fDensity;
}
- corsika::units::si::TimeType ArclengthFromGrammage(
+ corsika::units::si::LengthType ArclengthFromGrammage(
corsika::geometry::Trajectory<corsika::geometry::Line> const& pTraj,
corsika::units::si::GrammageType pGrammage) const override {
return pGrammage / fDensity;
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index b9e03eaa232a4132a077f585602100262b9c92d5..e1b44a3dc308720e7c370cc1500ab7056cc9d172 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -12,11 +12,12 @@
#ifndef _include_corsika_cascade_Cascade_h_
#define _include_corsika_cascade_Cascade_h_
+#include <corsika/environment/Environment.h>
#include <corsika/process/ProcessReturn.h>
#include <corsika/random/RNGManager.h>
+#include <corsika/random/UniformRealDistribution.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/units/PhysicalUnits.h>
-#include <corsika/environment/Environment.h>
#include <type_traits>
@@ -29,8 +30,10 @@ namespace corsika::cascade {
Cascade() = delete;
public:
- Cascade(Tracking& tr, ProcessList& pl, Stack& stack)
- : fTracking(tr)
+ Cascade(corsika::environment::Environment const& env, Tracking& tr, ProcessList& pl,
+ Stack& stack)
+ : fEnvironment(env)
+ , fTracking(tr)
, fProcessSequence(pl)
, fStack(stack) {}
@@ -53,26 +56,34 @@ namespace corsika::cascade {
}
void Step(Particle& particle) {
-
+ using namespace corsika::units::si;
// determine geometric tracking
corsika::setup::Trajectory step = fTracking.GetTrack(particle);
// determine combined total interaction length (inverse)
- InverseLengthType const total_inv_lambda =
+ InverseGrammageType const total_inv_lambda =
fProcessSequence.GetTotalInverseInteractionLength(particle, step);
- // sample random exponential step length
- std::exponential_distribution expDist(total_inv_lambda * 1_m);
- LengthType const next_interact = 1_m * expDist(fRNG);
+ // sample random exponential step length in grammage
+ std::exponential_distribution expDist((1_m * 1_m / 1_g) / total_inv_lambda);
+ GrammageType const next_interact = (1_g / (1_m * 1_m)) * expDist(fRNG);
std::cout << "total_inv_lambda=" << total_inv_lambda
<< ", next_interact=" << next_interact << std::endl;
+ // convert next_step from grammage to length
+ LengthType const distance_interact =
+ fEnvironment.GetUniverse()
+ ->GetContainingNode(particle.GetPosition())
+ ->GetModelProperties()
+ .ArclengthFromGrammage(step, next_interact);
+
// determine the maximum geometric step length
LengthType const distance_max = fProcessSequence.MaxStepLength(particle, step);
std::cout << "distance_max=" << distance_max << std::endl;
// determine combined total inverse decay time
- InverseTimeType const total_inv_lifetime = fProcessSequence.GetTotalInverseLifetime(particle);
+ InverseTimeType const total_inv_lifetime =
+ fProcessSequence.GetTotalInverseLifetime(particle);
// sample random exponential decay time
std::exponential_distribution expDistDecay(total_inv_lifetime * 1_s);
@@ -80,23 +91,20 @@ namespace corsika::cascade {
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 GrammageType distance_interact = fEnvironment.GetUniverse()->GetContainingNode(particle.GetPosition())->GetModelProperties().FromGrammage()
- // ....
-
// convert next_decay from time to length [m]
// Environment::GetDistance(step, next_decay);
- const double distance_decay = next_decay;
- // ....
+ LengthType const distance_decay = next_decay * particle.GetMomentum().norm() /
+ particle.GetEnergy() *
+ corsika::units::si::constants::cSquared;
// 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);
+ auto const min_distance =
+ std::min({distance_interact, distance_decay, 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));
+ particle.SetPosition(step.PositionFromArclength(min_distance));
+
// .... also update time, momentum, direction, ...
// apply all continuous processes on particle + track
@@ -107,27 +115,35 @@ namespace corsika::cascade {
// fStack.Delete(particle); // TODO: check if this is really needed
} else {
- 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) {
+ std::cout << "sth. happening before geometric limit ?"
+ << ((min_distance < distance_max) ? "yes" : "no") << std::endl;
+
+ if (min_distance < distance_max) { // interaction to happen within geometric limit
// check weather decay or interaction limits this step
- if (distance_decay > distance_interact) {
+
+ if (min_distance == distance_interact) {
std::cout << "collide" << std::endl;
- const double actual_inv_length =
+
+ InverseGrammageType const actual_inv_length =
fProcessSequence.GetTotalInverseInteractionLength(particle, step);
- const double sample_process = rmng() / (double)rmng.max();
- double inv_lambda_count = 0;
- fProcessSequence.SelectInteraction(particle, fStack, actual_inv_length,
- sample_process, inv_lambda_count);
+
+ corsika::random::UniformRealDistribution<InverseGrammageType> uniDist(
+ actual_inv_length);
+ const auto sample_process = uniDist(fRNG);
+ InverseGrammageType inv_lambda_count = 0. * meter * meter / gram;
+ fProcessSequence.SelectInteraction(particle, fStack, sample_process,
+ inv_lambda_count);
} else {
std::cout << "decay" << std::endl;
- const double actual_decay_time =
+ InverseTimeType const actual_decay_time =
fProcessSequence.GetTotalInverseLifetime(particle);
- const double sample_process = rmng() / (double)rmng.max();
- double inv_decay_count = 0;
- fProcessSequence.SelectDecay(particle, fStack, actual_decay_time, sample_process,
- inv_decay_count);
+
+ corsika::random::UniformRealDistribution<InverseTimeType> uniDist(
+ actual_decay_time);
+ const auto sample_process = uniDist(fRNG);
+ InverseTimeType inv_decay_count = 0 / second;
+ fProcessSequence.SelectDecay(particle, fStack, sample_process,
+ inv_decay_count);
}
}
}
@@ -139,8 +155,7 @@ namespace corsika::cascade {
Stack& fStack;
corsika::environment::Environment const& fEnvironment;
corsika::random::RNG& fRNG =
- corsika::random::RNGManager::GetInstance().GetRandomStream("cascade");
-
+ corsika::random::RNGManager::GetInstance().GetRandomStream("cascade");
};
} // namespace corsika::cascade
diff --git a/Framework/Cascade/testCascade.cc b/Framework/Cascade/testCascade.cc
index b39c3f145c697500141ad911bb110b8b7aba8da4..0e95d1dbdaea9573c4fc52fd25c9d667e31a8336 100644
--- a/Framework/Cascade/testCascade.cc
+++ b/Framework/Cascade/testCascade.cc
@@ -26,6 +26,10 @@
#include <corsika/geometry/RootCoordinateSystem.h>
#include <corsika/geometry/Vector.h>
+#include <corsika/environment/Environment.h>
+#include <corsika/environment/HomogeneousMedium.h>
+#include <corsika/environment/NuclearComposition.h>
+
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
using corsika::setup::Trajectory;
@@ -43,6 +47,27 @@ using namespace corsika::geometry;
using namespace std;
using namespace corsika::units::si;
+corsika::environment::Environment MakeDummyEnv() {
+ 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));
+
+ return env;
+}
+
static int fCount = 0;
class ProcessSplit : public corsika::process::ContinuousProcess<ProcessSplit> {
@@ -50,8 +75,8 @@ public:
ProcessSplit() {}
template <typename Particle, typename T>
- double MaxStepLength(Particle&, T&) const {
- return 1;
+ LengthType MaxStepLength(Particle&, T&) const {
+ return 1_m;
}
template <typename Particle, typename T, typename Stack>
@@ -81,16 +106,9 @@ private:
TEST_CASE("Cascade", "[Cascade]") {
corsika::random::RNGManager& rmng = corsika::random::RNGManager::GetInstance();
- rmng.RegisterRandomStream("s_rndm");
-
- 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));
+ rmng.RegisterRandomStream("cascade");
+ auto env = MakeDummyEnv();
tracking_line::TrackingLine<setup::Stack> tracking(env);
stack_inspector::StackInspector<setup::Stack> p0(true);
@@ -98,7 +116,7 @@ TEST_CASE("Cascade", "[Cascade]") {
const auto sequence = p0 + p1;
setup::Stack stack;
- corsika::cascade::Cascade EAS(tracking, sequence, stack);
+ corsika::cascade::Cascade EAS(env, tracking, sequence, stack);
CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
stack.Clear();
diff --git a/Framework/Geometry/Line.h b/Framework/Geometry/Line.h
index ae43908cf82cba819cb6110ad85f259370bbe98a..1b45f535a9c01acaa0084a42ca3b532360fa9b59 100644
--- a/Framework/Geometry/Line.h
+++ b/Framework/Geometry/Line.h
@@ -31,8 +31,10 @@ namespace corsika::geometry {
, v0(pV0) {}
Point GetPosition(corsika::units::si::TimeType t) const { return r0 + v0 * t; }
-
- Point PositionFromArclength(corsika::units::si::LengthType l) const { return r0 + v0.normalized() * l; }
+
+ Point PositionFromArclength(corsika::units::si::LengthType l) const {
+ return r0 + v0.normalized() * l;
+ }
LengthType ArcLength(corsika::units::si::TimeType t1,
corsika::units::si::TimeType t2) const {
diff --git a/Framework/Geometry/Trajectory.h b/Framework/Geometry/Trajectory.h
index 5bfcea0216ee6fd496275003b307cccc3f804ee8..0d142d91251ee15d22cb78c3daa8aeeb148b7699 100644
--- a/Framework/Geometry/Trajectory.h
+++ b/Framework/Geometry/Trajectory.h
@@ -15,9 +15,6 @@
#include <corsika/geometry/Point.h>
#include <corsika/units/PhysicalUnits.h>
-using corsika::units::si::LengthType;
-using corsika::units::si::TimeType;
-
namespace corsika::geometry {
template <typename T>
@@ -39,7 +36,7 @@ namespace corsika::geometry {
Point GetPosition(double u) const { return T::GetPosition(fTimeLength * u); }
- TimeType GetDuration() const { return fTimeLength; }
+ corsika::units::si::TimeType GetDuration() const { return fTimeLength; }
LengthType GetDistance(corsika::units::si::TimeType t) const {
assert(t > fTimeLength);
diff --git a/Framework/Geometry/testGeometry.cc b/Framework/Geometry/testGeometry.cc
index c51895f48f298f394c129c4498cbbe6db1179eed..2f1f891814cc68d7e5f585c9992dce63f244cb8a 100644
--- a/Framework/Geometry/testGeometry.cc
+++ b/Framework/Geometry/testGeometry.cc
@@ -164,11 +164,10 @@ TEST_CASE("Trajectories") {
QuantityVector<length_d>(4_m, 0_m, 0_m))
.norm()
.magnitude() == Approx(0).margin(absMargin));
-
- CHECK(
- (line.GetPosition(7_s) - line.PositionFromArclength(line.ArcLength(0_s, 7_s)))
- .norm()
- .magnitude() == Approx(0).margin(absMargin));
+
+ CHECK((line.GetPosition(7_s) - line.PositionFromArclength(line.ArcLength(0_s, 7_s)))
+ .norm()
+ .magnitude() == Approx(0).margin(absMargin));
auto const t = 1_s;
Trajectory<Line> base(line, t);
diff --git a/Framework/ProcessSequence/DecayProcess.h b/Framework/ProcessSequence/DecayProcess.h
index 8a49803ed49913a1c8cf5c217a5c866de307f4bd..076a32410495f8e7524b265284f7aa1d4f8417d5 100644
--- a/Framework/ProcessSequence/DecayProcess.h
+++ b/Framework/ProcessSequence/DecayProcess.h
@@ -14,6 +14,7 @@
#include <corsika/process/ProcessReturn.h> // for convenience
#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/units/PhysicalUnits.h>
namespace corsika::process {
@@ -38,10 +39,10 @@ namespace corsika::process {
inline EProcessReturn DoDecay(Particle&, Stack&) const;
template <typename Particle>
- inline double GetLifetime(Particle& p) const;
+ corsika::units::si::TimeType GetLifetime(Particle& p) const;
template <typename Particle>
- inline double GetInverseLifetime(Particle& p) const {
+ corsika::units::si::InverseTimeType GetInverseLifetime(Particle& p) const {
return 1. / GetRef().GetLifetime(p);
}
};
diff --git a/Framework/ProcessSequence/InteractionProcess.h b/Framework/ProcessSequence/InteractionProcess.h
index b0df1eebb5cc88cdb6a6b207906e0a1f2699bed1..a56edffa5f3def9dc866abead57142ec26db7890 100644
--- a/Framework/ProcessSequence/InteractionProcess.h
+++ b/Framework/ProcessSequence/InteractionProcess.h
@@ -14,6 +14,7 @@
#include <corsika/process/ProcessReturn.h> // for convenience
#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/units/PhysicalUnits.h>
namespace corsika::process {
@@ -37,11 +38,12 @@ namespace corsika::process {
template <typename P, typename S>
inline EProcessReturn DoInteraction(P&, S&) const;
- template <typename P, typename T>
- inline double GetInteractionLength(P&, T&) const;
+ template <typename Particle, typename Track>
+ corsika::units::si::GrammageType GetInteractionLength(Particle& p, Track& t) const;
template <typename Particle, typename Track>
- inline double GetInverseInteractionLength(Particle& p, Track& t) const {
+ corsika::units::si::InverseGrammageType GetInverseInteractionLength(Particle& p,
+ Track& t) const {
return 1. / GetRef().GetInteractionLength(p, t);
}
};
diff --git a/Framework/ProcessSequence/ProcessSequence.h b/Framework/ProcessSequence/ProcessSequence.h
index 3428b497eb10989689af94eb8ef9970ea3b82fbd..acf8a02f1c557d6b9977a5811c266c144d940908 100644
--- a/Framework/ProcessSequence/ProcessSequence.h
+++ b/Framework/ProcessSequence/ProcessSequence.h
@@ -18,6 +18,7 @@
//#include <corsika/process/DiscreteProcess.h>
#include <corsika/process/InteractionProcess.h>
#include <corsika/process/ProcessReturn.h>
+#include <corsika/units/PhysicalUnits.h>
#include <cmath>
#include <limits>
@@ -174,32 +175,41 @@ namespace corsika::process {
}
template <typename Particle, typename Track>
- LengthType MaxStepLength(Particle& p, Track& track) const {
- LengthType max_length = std::numeric_limits<double>::infinity() * 1_m;
+ corsika::units::si::LengthType MaxStepLength(Particle& p, Track& track) const {
+ corsika::units::si::LengthType max_length =
+ std::numeric_limits<double>::infinity() * corsika::units::si::meter;
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));
+ corsika::units::si::LengthType const len = A.MaxStepLength(p, track);
+ max_length = std::min(max_length, len);
}
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));
+ corsika::units::si::LengthType const len = B.MaxStepLength(p, track);
+ max_length = std::min(max_length, len);
}
return max_length;
}
template <typename Particle, typename Track>
- GrammageType GetTotalInteractionLength(Particle& p, Track& t) const {
+ corsika::units::si::GrammageType GetTotalInteractionLength(Particle& p,
+ Track& t) const {
return 1. / GetInverseInteractionLength(p, t);
}
template <typename Particle, typename Track>
- InverseGrammageType GetTotalInverseInteractionLength(Particle& p, Track& t) const {
+ corsika::units::si::InverseGrammageType GetTotalInverseInteractionLength(
+ Particle& p, Track& t) const {
return GetInverseInteractionLength(p, t);
}
template <typename Particle, typename Track>
- InverseGrammageType GetInverseInteractionLength(Particle& p, Track& t) const {
- InverseGrammageType tot = 0;
+ corsika::units::si::InverseGrammageType GetInverseInteractionLength(Particle& p,
+ Track& t) const {
+ using namespace corsika::units::si;
+
+ InverseGrammageType tot = 0 * meter * meter / gram;
+
if constexpr (std::is_base_of<InteractionProcess<T1>, T1>::value ||
is_process_sequence<T1>::value) {
tot += A.GetInverseInteractionLength(p, t);
@@ -212,21 +222,20 @@ namespace corsika::process {
}
template <typename Particle, typename Stack>
- inline EProcessReturn SelectInteraction(Particle& p, Stack& s,
- InverseGrammageType lambda_inv_tot,
- InverseGrammageType rndm_select,
- InverseGrammageType& lambda_inv_count) const {
+ EProcessReturn SelectInteraction(
+ Particle& p, Stack& s, corsika::units::si::InverseGrammageType lambda_select,
+ corsika::units::si::InverseGrammageType& 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);
+ A.SelectInteraction(p, s, 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<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_tot < lambda_inv_count) { // more pedagogical: rndm_select < lambda_inv_count / lambda_inv_tot
+ if (lambda_select < lambda_inv_count) {
A.DoInteraction(p, s);
return EProcessReturn::eInteracted;
}
@@ -235,14 +244,14 @@ namespace corsika::process {
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);
+ B.SelectInteraction(p, s, 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<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) {
+ if (lambda_select < lambda_inv_count) {
B.DoInteraction(p, s);
return EProcessReturn::eInteracted;
}
@@ -251,18 +260,21 @@ namespace corsika::process {
}
template <typename Particle>
- TimeType GetTotalLifetime(Particle& p) const {
+ corsika::units::si::TimeType GetTotalLifetime(Particle& p) const {
return 1. / GetInverseLifetime(p);
}
template <typename Particle>
- InverseTimeType GetTotalInverseLifetime(Particle& p) const {
+ corsika::units::si::InverseTimeType GetTotalInverseLifetime(Particle& p) const {
return GetInverseLifetime(p);
}
template <typename Particle>
- InverseTimeType GetInverseLifetime(Particle& p) const {
- InverseTimeType tot = 0;
+ corsika::units::si::InverseTimeType GetInverseLifetime(Particle& p) const {
+ using namespace corsika::units::si;
+
+ corsika::units::si::InverseTimeType tot = 0 / second;
+
if constexpr (std::is_base_of<DecayProcess<T1>, T1>::value ||
is_process_sequence<T1>::value) {
tot += A.GetInverseLifetime(p);
@@ -276,20 +288,20 @@ namespace corsika::process {
// select decay process
template <typename Particle, typename Stack>
- EProcessReturn SelectDecay(Particle& p, Stack& s, InverseTimeType decay_inv_tot,
- InverseTimeType rndm_select,
- InverseTimeType& decay_inv_count) const {
+ EProcessReturn SelectDecay(
+ Particle& p, Stack& s, corsika::units::si::InverseTimeType decay_select,
+ corsika::units::si::InverseTimeType& 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);
+ const EProcessReturn ret = A.SelectDecay(p, s, 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<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_tot < decay_inv_count) { // more pedagogical: rndm_select < decay_inv_count / decay_inv_tot
+ if (decay_select < decay_inv_count) { // more pedagogical: rndm_select <
+ // decay_inv_count / decay_inv_tot
A.DoDecay(p, s);
return EProcessReturn::eDecayed;
}
@@ -297,15 +309,14 @@ namespace corsika::process {
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);
+ const EProcessReturn ret = B.SelectDecay(p, s, 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<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) {
+ if (decay_select < decay_inv_count) {
B.DoDecay(p, s);
return EProcessReturn::eDecayed;
}
diff --git a/Framework/ProcessSequence/testProcessSequence.cc b/Framework/ProcessSequence/testProcessSequence.cc
index af9bfeeade6c2d149f88b803f0ddd2ebb105e6aa..e8a4828cda371a3bb67e0bd1de6b76b889292622 100644
--- a/Framework/ProcessSequence/testProcessSequence.cc
+++ b/Framework/ProcessSequence/testProcessSequence.cc
@@ -101,9 +101,9 @@ public:
return EProcessReturn::eOk;
}
template <typename Particle, typename Track>
- inline double GetInteractionLength(Particle&, Track&) const {
+ GrammageType GetInteractionLength(Particle&, Track&) const {
cout << "Process2::GetInteractionLength" << endl;
- return 3;
+ return 3_g / (1_cm * 1_cm);
}
};
@@ -124,9 +124,9 @@ public:
return EProcessReturn::eOk;
}
template <typename Particle, typename Track>
- inline double GetInteractionLength(Particle&, Track&) const {
+ GrammageType GetInteractionLength(Particle&, Track&) const {
cout << "Process3::GetInteractionLength" << endl;
- return 1.;
+ return 1_g / (1_cm * 1_cm);
}
};
@@ -165,8 +165,8 @@ public:
globalCount++;
}
template <typename Particle>
- double GetLifetime(Particle&) const {
- return 1;
+ TimeType GetLifetime(Particle&) const {
+ return 1_s;
}
template <typename Particle, typename Stack>
EProcessReturn DoDecay(Particle&, Stack&) const {
@@ -209,9 +209,9 @@ TEST_CASE("Process Sequence", "[Process Sequence]") {
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;
+ GrammageType const tot = sequence2.GetTotalInteractionLength(s, t);
+ InverseGrammageType const tot_inv = sequence2.GetTotalInverseInteractionLength(s, t);
+ cout << "lambda_tot=" << tot << "; lambda_tot_inv=" << tot_inv << endl;
}
SECTION("lifetime") {
@@ -223,9 +223,9 @@ TEST_CASE("Process Sequence", "[Process Sequence]") {
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;
+ TimeType const tot = sequence2.GetTotalLifetime(s);
+ InverseTimeType const tot_inv = sequence2.GetTotalInverseLifetime(s);
+ cout << "lambda_tot=" << tot << "; lambda_tot_inv=" << tot_inv << endl;
}
SECTION("sectionTwo") {
diff --git a/Framework/Units/PhysicalUnits.h b/Framework/Units/PhysicalUnits.h
index 734fd6f6ed76a6c15342ff0674b9d1dd996e3807..95af0e3d264e35a48a5e0e79d9370a619ea78e4e 100644
--- a/Framework/Units/PhysicalUnits.h
+++ b/Framework/Units/PhysicalUnits.h
@@ -54,9 +54,12 @@ namespace corsika::units::si {
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<area_d, double>;
- using InverseLengthType = phys::units::quantity<phys::units::dimensions<-1, 0, 0>, double>;
- using InverseTimeType = phys::units::quantity<phys::units::dimensions<0, 0, -1>, double>;
- using InverseGrammageType = phys::units::quantity<phys::units::dimensions<2, -1, 0>, double>;
+ using InverseLengthType =
+ phys::units::quantity<phys::units::dimensions<-1, 0, 0>, double>;
+ using InverseTimeType =
+ phys::units::quantity<phys::units::dimensions<0, 0, -1>, double>;
+ using InverseGrammageType =
+ phys::units::quantity<phys::units::dimensions<2, -1, 0>, double>;
} // end namespace corsika::units::si
@@ -75,7 +78,7 @@ namespace phys {
QUANTITY_DEFINE_SCALING_LITERALS(eV, energy_d,
magnitude(corsika::units::si::constants::eV))
- QUANTITY_DEFINE_SCALING_LITERALS(barn, corsika::units::si::sigma_d,
+ QUANTITY_DEFINE_SCALING_LITERALS(barn, corsika::units::si::area_d,
magnitude(corsika::units::si::constants::barn))
QUANTITY_DEFINE_SCALING_LITERALS(Ns, corsika::units::si::momentum_d,
diff --git a/Framework/Units/testUnits.cc b/Framework/Units/testUnits.cc
index 4990cd62462af914744aefe76be86528333a978e..db210a7738a57f373286904624766429e384cc56 100644
--- a/Framework/Units/testUnits.cc
+++ b/Framework/Units/testUnits.cc
@@ -43,8 +43,8 @@ TEST_CASE("PhysicalUnits", "[Units]") {
[[maybe_unused]] std::array<EnergyType, 4> arr3 = {1_GeV, 1_eV, 5_MeV};
- [[maybe_unused]] auto p1 = 10_newton_second;
- REQUIRE(p1 == 10_newton_second);
+ auto p1 = 10_Ns;
+ REQUIRE(p1 == 10_Ns);
}
SECTION("Powers in literal units") {
@@ -69,6 +69,8 @@ TEST_CASE("PhysicalUnits", "[Units]") {
REQUIRE(1_A / 1_EA == Approx(1e-18));
REQUIRE(1_K / 1_ZK == Approx(1e-21));
REQUIRE(1_mol / 1_Ymol == Approx(1e-24));
+
+ REQUIRE(std::min(1_A, 2_A) == 1_A);
}
SECTION("Powers and units") {
@@ -83,7 +85,7 @@ TEST_CASE("PhysicalUnits", "[Units]") {
const MassType m = 1_kg;
const SpeedType v = 1_m / 1_s;
- REQUIRE(m * v == 1_newton_second);
+ REQUIRE(m * v == 1_Ns);
const double lgE = log10(E2 / 1_GeV);
REQUIRE(lgE == Approx(log10(40.)));
diff --git a/Processes/StackInspector/StackInspector.cc b/Processes/StackInspector/StackInspector.cc
index ad90442848976c4cc779977a2c2f1890cedb1091..0dce10b08f9971c3265845a858fd8b66d943e831 100644
--- a/Processes/StackInspector/StackInspector.cc
+++ b/Processes/StackInspector/StackInspector.cc
@@ -58,8 +58,9 @@ process::EProcessReturn StackInspector<Stack>::DoContinuous(Particle&, setup::Tr
}
template <typename Stack>
-double StackInspector<Stack>::MaxStepLength(Particle&, setup::Trajectory&) const {
- return std::numeric_limits<double>::infinity();
+corsika::units::si::LengthType StackInspector<Stack>::MaxStepLength(
+ Particle&, setup::Trajectory&) const {
+ return std::numeric_limits<double>::infinity() * meter;
}
template <typename Stack>
diff --git a/Processes/StackInspector/StackInspector.h b/Processes/StackInspector/StackInspector.h
index 6032beac618b0941b01ffd76962c3bc13ad2300b..6a7d2fc62f248813895b3bfef4e673e010b20ad1 100644
--- a/Processes/StackInspector/StackInspector.h
+++ b/Processes/StackInspector/StackInspector.h
@@ -13,8 +13,8 @@
#define _Physics_StackInspector_StackInspector_h_
#include <corsika/process/ContinuousProcess.h>
-
#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/units/PhysicalUnits.h>
namespace corsika::process {
@@ -36,7 +36,8 @@ namespace corsika::process {
EProcessReturn DoContinuous(Particle&, corsika::setup::Trajectory&, Stack& s) const;
// template <typename Particle>
- double MaxStepLength(Particle&, corsika::setup::Trajectory&) const;
+ corsika::units::si::LengthType MaxStepLength(Particle&,
+ corsika::setup::Trajectory&) const;
private:
bool fReport;
diff --git a/Processes/TrackingLine/TrackingLine.h b/Processes/TrackingLine/TrackingLine.h
index 398d4c4e5e8ef389c6107076ae4ea2885bdd119d..fab03f0f0fc0993b79796e37c476ec0f01f8a423 100644
--- a/Processes/TrackingLine/TrackingLine.h
+++ b/Processes/TrackingLine/TrackingLine.h
@@ -76,12 +76,12 @@ namespace corsika::process {
auto GetTrack(Particle const& 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();
-
+
// to do: include effect of magnetic field
geometry::Line line(currentPosition, velocity);