diff --git a/Environment/HomogeneousMedium.h b/Environment/HomogeneousMedium.h
index a8fa5f7283cfe0e776383f24915b98748d54349f..019774d5120a605fade3c50c62c7baa88bb8eb00 100644
--- a/Environment/HomogeneousMedium.h
+++ b/Environment/HomogeneousMedium.h
@@ -43,15 +43,15 @@ namespace corsika::environment {
corsika::units::si::GrammageType IntegratedGrammage(
corsika::geometry::Trajectory<corsika::geometry::Line> const& pTraj,
- corsika::units::si::TimeType pTo) const override {
+ corsika::units::si::LengthType pTo) const override {
using namespace corsika::units::si;
- return pTraj.ArcLength(0_s, pTo) * fDensity;
+ pTo * fDensity;
}
- corsika::units::si::TimeType FromGrammage(
+ corsika::units::si::TimeType ArclengthFromGrammage(
corsika::geometry::Trajectory<corsika::geometry::Line> const& pTraj,
corsika::units::si::GrammageType pGrammage) const override {
- return pTraj.TimeFromArclength(pGrammage / fDensity);
+ return pGrammage / fDensity;
}
};
diff --git a/Environment/IMediumModel.h b/Environment/IMediumModel.h
index 4a1eedd86562187c310d180a4516ee3425a05674..d527c80435cb9ade5acfaf1c2bc9623fd44a1161 100644
--- a/Environment/IMediumModel.h
+++ b/Environment/IMediumModel.h
@@ -20,9 +20,9 @@ namespace corsika::environment {
// approach for now, only lines are supported
virtual corsika::units::si::GrammageType IntegratedGrammage(
corsika::geometry::Trajectory<corsika::geometry::Line> const&,
- corsika::units::si::TimeType) const = 0;
+ corsika::units::si::LengthType) const = 0;
- virtual corsika::units::si::TimeType FromGrammage(
+ virtual corsika::units::si::LengthType ArclengthFromGrammage(
corsika::geometry::Trajectory<corsika::geometry::Line> const&,
corsika::units::si::GrammageType) const = 0;
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index c311a5492541b419372111348a7c5be977427c29..b9e03eaa232a4132a077f585602100262b9c92d5 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -16,6 +16,7 @@
#include <corsika/random/RNGManager.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/units/PhysicalUnits.h>
+#include <corsika/environment/Environment.h>
#include <type_traits>
@@ -30,12 +31,12 @@ namespace corsika::cascade {
public:
Cascade(Tracking& tr, ProcessList& pl, Stack& stack)
: fTracking(tr)
- , fProcesseList(pl)
+ , fProcessSequence(pl)
, fStack(stack) {}
void Init() {
fTracking.Init();
- fProcesseList.Init();
+ fProcessSequence.Init();
fStack.Init();
}
@@ -58,7 +59,7 @@ namespace corsika::cascade {
// determine combined total interaction length (inverse)
InverseLengthType const total_inv_lambda =
- fProcesseList.GetTotalInverseInteractionLength(particle, step);
+ fProcessSequence.GetTotalInverseInteractionLength(particle, step);
// sample random exponential step length
std::exponential_distribution expDist(total_inv_lambda * 1_m);
@@ -67,11 +68,11 @@ namespace corsika::cascade {
<< ", next_interact=" << next_interact << std::endl;
// determine the maximum geometric step length
- LengthType const distance_max = fProcesseList.MaxStepLength(particle, step);
+ 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 = fProcesseList.GetTotalInverseLifetime(particle);
+ InverseTimeType const total_inv_lifetime = fProcessSequence.GetTotalInverseLifetime(particle);
// sample random exponential decay time
std::exponential_distribution expDistDecay(total_inv_lifetime * 1_s);
@@ -81,7 +82,7 @@ namespace corsika::cascade {
// convert next_step from grammage to length [m]
// Environment::GetDistance(step, next_step);
- const GrammageType distance_interact = ;
+ const GrammageType distance_interact = fEnvironment.GetUniverse()->GetContainingNode(particle.GetPosition())->GetModelProperties().FromGrammage()
// ....
// convert next_decay from time to length [m]
@@ -93,14 +94,14 @@ namespace corsika::cascade {
const double distance_decay_interact = std::min(next_decay, next_interact);
const double distance_next = std::min(distance_decay_interact, distance_max);
- /// here the particle is actually moved along the trajectory to new position:
+ // here the particle is actually moved along the trajectory to new position:
// std::visit(corsika::setup::ParticleUpdate<Particle>{particle}, step);
particle.SetPosition(step.GetPosition(1));
// .... also update time, momentum, direction, ...
// apply all continuous processes on particle + track
corsika::process::EProcessReturn status =
- fProcesseList.DoContinuous(particle, step, fStack);
+ fProcessSequence.DoContinuous(particle, step, fStack);
if (status == corsika::process::EProcessReturn::eParticleAbsorbed) {
// fStack.Delete(particle); // TODO: check if this is really needed
@@ -114,18 +115,18 @@ namespace corsika::cascade {
if (distance_decay > distance_interact) {
std::cout << "collide" << std::endl;
const double actual_inv_length =
- fProcesseList.GetTotalInverseInteractionLength(particle, step);
+ fProcessSequence.GetTotalInverseInteractionLength(particle, step);
const double sample_process = rmng() / (double)rmng.max();
double inv_lambda_count = 0;
- fProcesseList.SelectInteraction(particle, fStack, actual_inv_length,
+ fProcessSequence.SelectInteraction(particle, fStack, actual_inv_length,
sample_process, inv_lambda_count);
} else {
std::cout << "decay" << std::endl;
const double actual_decay_time =
- fProcesseList.GetTotalInverseLifetime(particle);
+ fProcessSequence.GetTotalInverseLifetime(particle);
const double sample_process = rmng() / (double)rmng.max();
double inv_decay_count = 0;
- fProcesseList.SelectDecay(particle, fStack, actual_decay_time, sample_process,
+ fProcessSequence.SelectDecay(particle, fStack, actual_decay_time, sample_process,
inv_decay_count);
}
}
@@ -134,7 +135,7 @@ namespace corsika::cascade {
private:
Tracking& fTracking;
- ProcessList& fProcesseList;
+ ProcessList& fProcessSequence;
Stack& fStack;
corsika::environment::Environment const& fEnvironment;
corsika::random::RNG& fRNG =
diff --git a/Framework/ProcessSequence/ProcessSequence.h b/Framework/ProcessSequence/ProcessSequence.h
index 8517d4345202f8fb667f830eab3c093937b2ad05..c9a23fdc218e229a6eb3940ae71d76ea8da9a584 100644
--- a/Framework/ProcessSequence/ProcessSequence.h
+++ b/Framework/ProcessSequence/ProcessSequence.h
@@ -160,7 +160,7 @@ namespace corsika::process {
// void Hello() const { detail::CallHello<T1,T2>::Call(A, B); }
template <typename Particle, typename Track, typename Stack>
- inline EProcessReturn DoContinuous(Particle& p, Track& t, Stack& s) const {
+ EProcessReturn DoContinuous(Particle& p, Track& t, Stack& s) const {
EProcessReturn ret = EProcessReturn::eOk;
if constexpr (std::is_base_of<ContinuousProcess<T1>, T1>::value ||
is_process_sequence<T1>::value) {
@@ -174,8 +174,8 @@ namespace corsika::process {
}
template <typename Particle, typename Track>
- inline double MaxStepLength(Particle& p, Track& track) const {
- double max_length = std::numeric_limits<double>::infinity();
+ LengthType MaxStepLength(Particle& p, Track& track) const {
+ LengthType max_length = std::numeric_limits<double>::infinity() * 1_m;
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));
@@ -188,18 +188,18 @@ namespace corsika::process {
}
template <typename Particle, typename Track>
- inline double GetTotalInteractionLength(Particle& p, Track& t) const {
+ GrammageType GetTotalInteractionLength(Particle& p, Track& t) const {
return 1. / GetInverseInteractionLength(p, t);
}
template <typename Particle, typename Track>
- inline double GetTotalInverseInteractionLength(Particle& p, Track& t) const {
+ InverseGrammageType GetTotalInverseInteractionLength(Particle& p, Track& t) const {
return GetInverseInteractionLength(p, t);
}
template <typename Particle, typename Track>
- inline double GetInverseInteractionLength(Particle& p, Track& t) const {
- double tot = 0;
+ InverseGrammageType GetInverseInteractionLength(Particle& p, Track& t) const {
+ InverseGrammageType tot = 0;
if constexpr (std::is_base_of<InteractionProcess<T1>, T1>::value ||
is_process_sequence<T1>::value) {
tot += A.GetInverseInteractionLength(p, t);
@@ -213,9 +213,9 @@ namespace corsika::process {
template <typename Particle, typename Stack>
inline EProcessReturn SelectInteraction(Particle& p, Stack& s,
- const double lambda_inv_tot,
- const double rndm_select,
- double& lambda_inv_count) const {
+ InverseGrammageType lambda_inv_tot,
+ InverseGrammageType rndm_select,
+ InverseGrammageType& lambda_inv_count) const {
if constexpr (is_process_sequence<T1>::value) {
// if A is a process sequence --> check inside
const EProcessReturn ret =
@@ -226,7 +226,7 @@ namespace corsika::process {
// if this is not a ContinuousProcess --> evaluate probability
lambda_inv_count += A.GetInverseInteractionLength(p, s);
// check if we should execute THIS process and then EXIT
- if (rndm_select < lambda_inv_count / lambda_inv_tot) {
+ if (rndm_select * lambda_inv_tot < lambda_inv_count) { // more pedagogical: rndm_select < lambda_inv_count / lambda_inv_tot
A.DoInteraction(p, s);
return EProcessReturn::eInteracted;
}
@@ -251,18 +251,18 @@ namespace corsika::process {
}
template <typename Particle>
- inline double GetTotalLifetime(Particle& p) const {
+ TimeType GetTotalLifetime(Particle& p) const {
return 1. / GetInverseLifetime(p);
}
template <typename Particle>
- inline double GetTotalInverseLifetime(Particle& p) const {
+ InverseTimeType GetTotalInverseLifetime(Particle& p) const {
return GetInverseLifetime(p);
}
template <typename Particle>
- inline double GetInverseLifetime(Particle& p) const {
- double tot = 0;
+ InverseTimeType GetInverseLifetime(Particle& p) const {
+ InverseTimeType tot = 0;
if constexpr (std::is_base_of<DecayProcess<T1>, T1>::value ||
is_process_sequence<T1>::value) {
tot += A.GetInverseLifetime(p);
@@ -276,9 +276,9 @@ namespace corsika::process {
// select decay process
template <typename Particle, typename Stack>
- inline EProcessReturn SelectDecay(Particle& p, Stack& s, const double decay_inv_tot,
- const double rndm_select,
- double& decay_inv_count) const {
+ EProcessReturn SelectDecay(Particle& p, Stack& s, InverseTimeType decay_inv_tot,
+ InverseTimeType rndm_select,
+ InverseTimeType& decay_inv_count) const {
if constexpr (is_process_sequence<T1>::value) {
// if A is a process sequence --> check inside
const EProcessReturn ret =
@@ -289,7 +289,7 @@ namespace corsika::process {
// if this is not a ContinuousProcess --> evaluate probability
decay_inv_count += A.GetInverseLifetime(p);
// check if we should execute THIS process and then EXIT
- if (rndm_select < decay_inv_count / decay_inv_tot) {
+ if (rndm_select * decay_inv_tot < decay_inv_count) { // more pedagogical: rndm_select < decay_inv_count / decay_inv_tot
A.DoDecay(p, s);
return EProcessReturn::eDecayed;
}
@@ -314,7 +314,7 @@ namespace corsika::process {
}
/// TODO the const_cast is not nice, think about the constness here
- inline void Init() const {
+ void Init() const {
const_cast<T1*>(&A)->Init();
const_cast<T2*>(&B)->Init();
}
diff --git a/Framework/Units/PhysicalUnits.h b/Framework/Units/PhysicalUnits.h
index 7ebe8fe0d792f6f3f4767996ae5684f392e28ca9..734fd6f6ed76a6c15342ff0674b9d1dd996e3807 100644
--- a/Framework/Units/PhysicalUnits.h
+++ b/Framework/Units/PhysicalUnits.h
@@ -56,6 +56,7 @@ namespace corsika::units::si {
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>;
} // end namespace corsika::units::si