IAP GITLAB

Skip to content
Snippets Groups Projects
Commit e9023467 authored by Maximilian Reininghaus's avatar Maximilian Reininghaus :vulcan:
Browse files

changed some units

parent aa7683f7
No related branches found
No related tags found
1 merge request!29Resolve "Consistent units in cascade step"
Pipeline #58 failed
Stage: build
Hide whitespace changes
Inline Side-by-side
Environment/HomogeneousMedium.h
+ 4
− 4
View file @ e9023467
...@@ -43,15 +43,15 @@ namespace corsika::environment { ...@@ -43,15 +43,15 @@ namespace corsika::environment {
corsika::units::si::GrammageType IntegratedGrammage( corsika::units::si::GrammageType IntegratedGrammage(
corsika::geometry::Trajectory<corsika::geometry::Line> const& pTraj, 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; 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::geometry::Trajectory<corsika::geometry::Line> const& pTraj,
corsika::units::si::GrammageType pGrammage) const override { corsika::units::si::GrammageType pGrammage) const override {
return pTraj.TimeFromArclength(pGrammage / fDensity); return pGrammage / fDensity;
} }
}; };
......
Environment/IMediumModel.h
+ 2
− 2
View file @ e9023467
...@@ -20,9 +20,9 @@ namespace corsika::environment { ...@@ -20,9 +20,9 @@ namespace corsika::environment {
// approach for now, only lines are supported // approach for now, only lines are supported
virtual corsika::units::si::GrammageType IntegratedGrammage( virtual corsika::units::si::GrammageType IntegratedGrammage(
corsika::geometry::Trajectory<corsika::geometry::Line> const&, 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::geometry::Trajectory<corsika::geometry::Line> const&,
corsika::units::si::GrammageType) const = 0; corsika::units::si::GrammageType) const = 0;
......
Framework/Cascade/Cascade.h
+ 14
− 13
View file @ e9023467
...@@ -16,6 +16,7 @@ ...@@ -16,6 +16,7 @@
#include <corsika/random/RNGManager.h> #include <corsika/random/RNGManager.h>
#include <corsika/setup/SetupTrajectory.h> #include <corsika/setup/SetupTrajectory.h>
#include <corsika/units/PhysicalUnits.h> #include <corsika/units/PhysicalUnits.h>
#include <corsika/environment/Environment.h>
#include <type_traits> #include <type_traits>
...@@ -30,12 +31,12 @@ namespace corsika::cascade { ...@@ -30,12 +31,12 @@ namespace corsika::cascade {
public: public:
Cascade(Tracking& tr, ProcessList& pl, Stack& stack) Cascade(Tracking& tr, ProcessList& pl, Stack& stack)
: fTracking(tr) : fTracking(tr)
, fProcesseList(pl) , fProcessSequence(pl)
, fStack(stack) {} , fStack(stack) {}
void Init() { void Init() {
fTracking.Init(); fTracking.Init();
fProcesseList.Init(); fProcessSequence.Init();
fStack.Init(); fStack.Init();
} }
...@@ -58,7 +59,7 @@ namespace corsika::cascade { ...@@ -58,7 +59,7 @@ namespace corsika::cascade {
// determine combined total interaction length (inverse) // determine combined total interaction length (inverse)
InverseLengthType const total_inv_lambda = InverseLengthType const total_inv_lambda =
fProcesseList.GetTotalInverseInteractionLength(particle, step); fProcessSequence.GetTotalInverseInteractionLength(particle, step);
// sample random exponential step length // sample random exponential step length
std::exponential_distribution expDist(total_inv_lambda * 1_m); std::exponential_distribution expDist(total_inv_lambda * 1_m);
...@@ -67,11 +68,11 @@ namespace corsika::cascade { ...@@ -67,11 +68,11 @@ namespace corsika::cascade {
<< ", next_interact=" << next_interact << std::endl; << ", next_interact=" << next_interact << std::endl;
// determine the maximum geometric step length // 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; std::cout << "distance_max=" << distance_max << std::endl;
// determine combined total inverse decay time // 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 // sample random exponential decay time
std::exponential_distribution expDistDecay(total_inv_lifetime * 1_s); std::exponential_distribution expDistDecay(total_inv_lifetime * 1_s);
...@@ -81,7 +82,7 @@ namespace corsika::cascade { ...@@ -81,7 +82,7 @@ namespace corsika::cascade {
// convert next_step from grammage to length [m] // convert next_step from grammage to length [m]
// Environment::GetDistance(step, next_step); // 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] // convert next_decay from time to length [m]
...@@ -93,14 +94,14 @@ namespace corsika::cascade { ...@@ -93,14 +94,14 @@ namespace corsika::cascade {
const double distance_decay_interact = std::min(next_decay, next_interact); const double distance_decay_interact = std::min(next_decay, next_interact);
const double distance_next = std::min(distance_decay_interact, distance_max); 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); // std::visit(corsika::setup::ParticleUpdate<Particle>{particle}, step);
particle.SetPosition(step.GetPosition(1)); particle.SetPosition(step.GetPosition(1));
// .... also update time, momentum, direction, ... // .... also update time, momentum, direction, ...
// apply all continuous processes on particle + track // apply all continuous processes on particle + track
corsika::process::EProcessReturn status = corsika::process::EProcessReturn status =
fProcesseList.DoContinuous(particle, step, fStack); fProcessSequence.DoContinuous(particle, step, fStack);
if (status == corsika::process::EProcessReturn::eParticleAbsorbed) { if (status == corsika::process::EProcessReturn::eParticleAbsorbed) {
// fStack.Delete(particle); // TODO: check if this is really needed // fStack.Delete(particle); // TODO: check if this is really needed
...@@ -114,18 +115,18 @@ namespace corsika::cascade { ...@@ -114,18 +115,18 @@ namespace corsika::cascade {
if (distance_decay > distance_interact) { if (distance_decay > distance_interact) {
std::cout << "collide" << std::endl; std::cout << "collide" << std::endl;
const double actual_inv_length = const double actual_inv_length =
fProcesseList.GetTotalInverseInteractionLength(particle, step); fProcessSequence.GetTotalInverseInteractionLength(particle, step);
const double sample_process = rmng() / (double)rmng.max(); const double sample_process = rmng() / (double)rmng.max();
double inv_lambda_count = 0; double inv_lambda_count = 0;
fProcesseList.SelectInteraction(particle, fStack, actual_inv_length, fProcessSequence.SelectInteraction(particle, fStack, actual_inv_length,
sample_process, inv_lambda_count); sample_process, inv_lambda_count);
} else { } else {
std::cout << "decay" << std::endl; std::cout << "decay" << std::endl;
const double actual_decay_time = const double actual_decay_time =
fProcesseList.GetTotalInverseLifetime(particle); fProcessSequence.GetTotalInverseLifetime(particle);
const double sample_process = rmng() / (double)rmng.max(); const double sample_process = rmng() / (double)rmng.max();
double inv_decay_count = 0; 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); inv_decay_count);
} }
} }
...@@ -134,7 +135,7 @@ namespace corsika::cascade { ...@@ -134,7 +135,7 @@ namespace corsika::cascade {
private: private:
Tracking& fTracking; Tracking& fTracking;
ProcessList& fProcesseList; ProcessList& fProcessSequence;
Stack& fStack; Stack& fStack;
corsika::environment::Environment const& fEnvironment; corsika::environment::Environment const& fEnvironment;
corsika::random::RNG& fRNG = corsika::random::RNG& fRNG =
......
Framework/ProcessSequence/ProcessSequence.h
+ 20
− 20
View file @ e9023467
...@@ -160,7 +160,7 @@ namespace corsika::process { ...@@ -160,7 +160,7 @@ namespace corsika::process {
// void Hello() const { detail::CallHello<T1,T2>::Call(A, B); } // void Hello() const { detail::CallHello<T1,T2>::Call(A, B); }
template <typename Particle, typename Track, typename Stack> 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; EProcessReturn ret = EProcessReturn::eOk;
if constexpr (std::is_base_of<ContinuousProcess<T1>, T1>::value || if constexpr (std::is_base_of<ContinuousProcess<T1>, T1>::value ||
is_process_sequence<T1>::value) { is_process_sequence<T1>::value) {
...@@ -174,8 +174,8 @@ namespace corsika::process { ...@@ -174,8 +174,8 @@ namespace corsika::process {
} }
template <typename Particle, typename Track> template <typename Particle, typename Track>
inline double MaxStepLength(Particle& p, Track& track) const { LengthType MaxStepLength(Particle& p, Track& track) const {
double max_length = std::numeric_limits<double>::infinity(); LengthType max_length = std::numeric_limits<double>::infinity() * 1_m;
if constexpr (std::is_base_of<ContinuousProcess<T1>, T1>::value || if constexpr (std::is_base_of<ContinuousProcess<T1>, T1>::value ||
is_process_sequence<T1>::value) { is_process_sequence<T1>::value) {
max_length = std::min(max_length, A.MaxStepLength(p, track)); max_length = std::min(max_length, A.MaxStepLength(p, track));
...@@ -188,18 +188,18 @@ namespace corsika::process { ...@@ -188,18 +188,18 @@ namespace corsika::process {
} }
template <typename Particle, typename Track> 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); return 1. / GetInverseInteractionLength(p, t);
} }
template <typename Particle, typename Track> template <typename Particle, typename Track>
inline double GetTotalInverseInteractionLength(Particle& p, Track& t) const { InverseGrammageType GetTotalInverseInteractionLength(Particle& p, Track& t) const {
return GetInverseInteractionLength(p, t); return GetInverseInteractionLength(p, t);
} }
template <typename Particle, typename Track> template <typename Particle, typename Track>
inline double GetInverseInteractionLength(Particle& p, Track& t) const { InverseGrammageType GetInverseInteractionLength(Particle& p, Track& t) const {
double tot = 0; InverseGrammageType tot = 0;
if constexpr (std::is_base_of<InteractionProcess<T1>, T1>::value || if constexpr (std::is_base_of<InteractionProcess<T1>, T1>::value ||
is_process_sequence<T1>::value) { is_process_sequence<T1>::value) {
tot += A.GetInverseInteractionLength(p, t); tot += A.GetInverseInteractionLength(p, t);
...@@ -213,9 +213,9 @@ namespace corsika::process { ...@@ -213,9 +213,9 @@ namespace corsika::process {
template <typename Particle, typename Stack> template <typename Particle, typename Stack>
inline EProcessReturn SelectInteraction(Particle& p, Stack& s, inline EProcessReturn SelectInteraction(Particle& p, Stack& s,
const double lambda_inv_tot, InverseGrammageType lambda_inv_tot,
const double rndm_select, InverseGrammageType rndm_select,
double& lambda_inv_count) const { InverseGrammageType& lambda_inv_count) const {
if constexpr (is_process_sequence<T1>::value) { if constexpr (is_process_sequence<T1>::value) {
// if A is a process sequence --> check inside // if A is a process sequence --> check inside
const EProcessReturn ret = const EProcessReturn ret =
...@@ -226,7 +226,7 @@ namespace corsika::process { ...@@ -226,7 +226,7 @@ namespace corsika::process {
// if this is not a ContinuousProcess --> evaluate probability // if this is not a ContinuousProcess --> evaluate probability
lambda_inv_count += A.GetInverseInteractionLength(p, s); lambda_inv_count += A.GetInverseInteractionLength(p, s);
// check if we should execute THIS process and then EXIT // 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); A.DoInteraction(p, s);
return EProcessReturn::eInteracted; return EProcessReturn::eInteracted;
} }
...@@ -251,18 +251,18 @@ namespace corsika::process { ...@@ -251,18 +251,18 @@ namespace corsika::process {
} }
template <typename Particle> template <typename Particle>
inline double GetTotalLifetime(Particle& p) const { TimeType GetTotalLifetime(Particle& p) const {
return 1. / GetInverseLifetime(p); return 1. / GetInverseLifetime(p);
} }
template <typename Particle> template <typename Particle>
inline double GetTotalInverseLifetime(Particle& p) const { InverseTimeType GetTotalInverseLifetime(Particle& p) const {
return GetInverseLifetime(p); return GetInverseLifetime(p);
} }
template <typename Particle> template <typename Particle>
inline double GetInverseLifetime(Particle& p) const { InverseTimeType GetInverseLifetime(Particle& p) const {
double tot = 0; InverseTimeType tot = 0;
if constexpr (std::is_base_of<DecayProcess<T1>, T1>::value || if constexpr (std::is_base_of<DecayProcess<T1>, T1>::value ||
is_process_sequence<T1>::value) { is_process_sequence<T1>::value) {
tot += A.GetInverseLifetime(p); tot += A.GetInverseLifetime(p);
...@@ -276,9 +276,9 @@ namespace corsika::process { ...@@ -276,9 +276,9 @@ namespace corsika::process {
// select decay process // select decay process
template <typename Particle, typename Stack> template <typename Particle, typename Stack>
inline EProcessReturn SelectDecay(Particle& p, Stack& s, const double decay_inv_tot, EProcessReturn SelectDecay(Particle& p, Stack& s, InverseTimeType decay_inv_tot,
const double rndm_select, InverseTimeType rndm_select,
double& decay_inv_count) const { InverseTimeType& decay_inv_count) const {
if constexpr (is_process_sequence<T1>::value) { if constexpr (is_process_sequence<T1>::value) {
// if A is a process sequence --> check inside // if A is a process sequence --> check inside
const EProcessReturn ret = const EProcessReturn ret =
...@@ -289,7 +289,7 @@ namespace corsika::process { ...@@ -289,7 +289,7 @@ namespace corsika::process {
// if this is not a ContinuousProcess --> evaluate probability // if this is not a ContinuousProcess --> evaluate probability
decay_inv_count += A.GetInverseLifetime(p); decay_inv_count += A.GetInverseLifetime(p);
// check if we should execute THIS process and then EXIT // 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); A.DoDecay(p, s);
return EProcessReturn::eDecayed; return EProcessReturn::eDecayed;
} }
...@@ -314,7 +314,7 @@ namespace corsika::process { ...@@ -314,7 +314,7 @@ namespace corsika::process {
} }
/// TODO the const_cast is not nice, think about the constness here /// 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<T1*>(&A)->Init();
const_cast<T2*>(&B)->Init(); const_cast<T2*>(&B)->Init();
} }
......
Framework/Units/PhysicalUnits.h
+ 1
− 0
View file @ e9023467
...@@ -56,6 +56,7 @@ namespace corsika::units::si { ...@@ -56,6 +56,7 @@ namespace corsika::units::si {
using CrossSectionType = phys::units::quantity<area_d, double>; using CrossSectionType = phys::units::quantity<area_d, double>;
using InverseLengthType = phys::units::quantity<phys::units::dimensions<-1, 0, 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 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 } // end namespace corsika::units::si
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

data privacy - impressum