From fb2adb40350f252acc3bb15b5e4be84d13dcde81 Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@kit.edu>
Date: Sun, 16 Dec 2018 16:27:54 +0100
Subject: [PATCH] towards integration of environment in cascade
---
Framework/Cascade/Cascade.h | 30 ++++++++++++---------------
Framework/Geometry/Helix.h | 8 +++++--
Framework/Geometry/Line.h | 2 ++
Framework/Geometry/testGeometry.cc | 26 ++++++++++++++++++-----
Framework/Units/PhysicalUnits.h | 14 ++++++-------
Processes/TrackingLine/TrackingLine.h | 11 ++++++----
6 files changed, 55 insertions(+), 36 deletions(-)
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index 7cd9d077d..c311a5492 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -19,15 +19,11 @@
#include <type_traits>
-using namespace corsika;
-using namespace corsika::units::si;
-
namespace corsika::cascade {
template <typename Tracking, typename ProcessList, typename Stack>
class Cascade {
-
- typedef typename Stack::ParticleType Particle;
+ using Particle = typename Stack::ParticleType;
Cascade() = delete;
@@ -57,39 +53,35 @@ namespace corsika::cascade {
void Step(Particle& particle) {
- // get access to random number generator
- static corsika::random::RNG& rmng =
- corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
-
// determine geometric tracking
corsika::setup::Trajectory step = fTracking.GetTrack(particle);
// determine combined total interaction length (inverse)
- const double total_inv_lambda =
+ InverseLengthType const total_inv_lambda =
fProcesseList.GetTotalInverseInteractionLength(particle, step);
// sample random exponential step length
- const double sample_step = rmng() / (double)rmng.max();
- const double next_interact = -log(sample_step) / total_inv_lambda;
+ std::exponential_distribution expDist(total_inv_lambda * 1_m);
+ LengthType const next_interact = 1_m * expDist(fRNG);
std::cout << "total_inv_lambda=" << total_inv_lambda
<< ", next_interact=" << next_interact << std::endl;
// determine the maximum geometric step length
- const double distance_max = fProcesseList.MaxStepLength(particle, step);
+ LengthType const distance_max = fProcesseList.MaxStepLength(particle, step);
std::cout << "distance_max=" << distance_max << std::endl;
// determine combined total inverse decay time
- const double total_inv_lifetime = fProcesseList.GetTotalInverseLifetime(particle);
+ InverseTimeType const total_inv_lifetime = fProcesseList.GetTotalInverseLifetime(particle);
// sample random exponential decay time
- const double sample_decay = rmng() / (double)rmng.max();
- const double next_decay = -log(sample_decay) / total_inv_lifetime;
+ std::exponential_distribution expDistDecay(total_inv_lifetime * 1_s);
+ TimeType const next_decay = 1_s * expDistDecay(fRNG);
std::cout << "total_inv_lifetime=" << total_inv_lifetime
<< ", next_decay=" << next_decay << std::endl;
// convert next_step from grammage to length [m]
// Environment::GetDistance(step, next_step);
- const double distance_interact = next_interact;
+ const GrammageType distance_interact = ;
// ....
// convert next_decay from time to length [m]
@@ -144,6 +136,10 @@ namespace corsika::cascade {
Tracking& fTracking;
ProcessList& fProcesseList;
Stack& fStack;
+ corsika::environment::Environment const& fEnvironment;
+ corsika::random::RNG& fRNG =
+ corsika::random::RNGManager::GetInstance().GetRandomStream("cascade");
+
};
} // namespace corsika::cascade
diff --git a/Framework/Geometry/Helix.h b/Framework/Geometry/Helix.h
index 2c8f6b879..27e3575d3 100644
--- a/Framework/Geometry/Helix.h
+++ b/Framework/Geometry/Helix.h
@@ -56,6 +56,10 @@ namespace corsika::geometry {
(vPerp * (cos(omegaC * t) - 1) + uPerp * sin(omegaC * t)) / omegaC;
}
+ Point PositionFromArclength(corsika::units::si::LengthType l) const {
+ return GetPosition(TimeFromArclength(l));
+ }
+
auto GetRadius() const { return radius; }
corsika::units::si::LengthType ArcLength(corsika::units::si::TimeType t1,
@@ -64,8 +68,8 @@ namespace corsika::geometry {
}
corsika::units::si::TimeType TimeFromArclength(
- corsika::units::si::LengthType t) const {
- return t / (vPar + vPerp).norm();
+ corsika::units::si::LengthType l) const {
+ return l / (vPar + vPerp).norm();
}
};
diff --git a/Framework/Geometry/Line.h b/Framework/Geometry/Line.h
index a79ba9ce8..ae43908cf 100644
--- a/Framework/Geometry/Line.h
+++ b/Framework/Geometry/Line.h
@@ -31,6 +31,8 @@ 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; }
LengthType ArcLength(corsika::units::si::TimeType t1,
corsika::units::si::TimeType t2) const {
diff --git a/Framework/Geometry/testGeometry.cc b/Framework/Geometry/testGeometry.cc
index 8130114aa..c51895f48 100644
--- a/Framework/Geometry/testGeometry.cc
+++ b/Framework/Geometry/testGeometry.cc
@@ -152,11 +152,21 @@ TEST_CASE("Trajectories") {
SECTION("Line") {
Vector<SpeedType::dimension_type> v0(rootCS,
- {1_m / second, 0_m / second, 0_m / second});
+ {3_m / second, 0_m / second, 0_m / second});
Line const line(r0, v0);
CHECK(
- (line.GetPosition(2_s).GetCoordinates() - QuantityVector<length_d>(2_m, 0_m, 0_m))
+ (line.GetPosition(2_s).GetCoordinates() - QuantityVector<length_d>(6_m, 0_m, 0_m))
+ .norm()
+ .magnitude() == Approx(0).margin(absMargin));
+
+ CHECK((line.PositionFromArclength(4_m).GetCoordinates() -
+ 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));
@@ -164,7 +174,7 @@ TEST_CASE("Trajectories") {
Trajectory<Line> base(line, t);
CHECK(line.GetPosition(t).GetCoordinates() == base.GetPosition(1.).GetCoordinates());
- CHECK(base.ArcLength(1_s, 2_s) / 1_m == Approx(1));
+ CHECK(base.ArcLength(1_s, 2_s) / 1_m == Approx(3));
}
SECTION("Helix") {
@@ -174,7 +184,8 @@ TEST_CASE("Trajectories") {
Vector<SpeedType::dimension_type> const vPerp(
rootCS, {3_m / second, 0_m / second, 0_m / second});
- auto const omegaC = 2 * M_PI / 1_s;
+ auto const T = 1_s;
+ auto const omegaC = 2 * M_PI / T;
Helix const helix(r0, omegaC, vPar, vPerp);
@@ -188,10 +199,15 @@ TEST_CASE("Trajectories") {
.norm()
.magnitude() == Approx(0).margin(absMargin));
+ CHECK(
+ (helix.GetPosition(7_s) - helix.PositionFromArclength(helix.ArcLength(0_s, 7_s)))
+ .norm()
+ .magnitude() == Approx(0).margin(absMargin));
+
auto const t = 1234_s;
Trajectory<Helix> const base(helix, t);
CHECK(helix.GetPosition(t).GetCoordinates() == base.GetPosition(1.).GetCoordinates());
- CHECK(base.ArcLength(1_s, 2_s) / 1_m == Approx(5));
+ CHECK(base.ArcLength(0_s, 1_s) / 1_m == Approx(5));
}
}
diff --git a/Framework/Units/PhysicalUnits.h b/Framework/Units/PhysicalUnits.h
index 7083ce743..7ebe8fe0d 100644
--- a/Framework/Units/PhysicalUnits.h
+++ b/Framework/Units/PhysicalUnits.h
@@ -39,8 +39,7 @@ namespace corsika::units::si {
constexpr phys::units::quantity<momentum_d> newton_second{meter * kilogram / second};
/// defining cross section
- using sigma_d = phys::units::dimensions<2, 0, 0>;
- constexpr phys::units::quantity<sigma_d> barn{Rep(1.e-28L) * meter * meter};
+ constexpr phys::units::quantity<area_d> barn{Rep(1.e-28L) * meter * meter};
/// add the unit-types
using LengthType = phys::units::quantity<phys::units::length_d, double>;
@@ -54,7 +53,9 @@ namespace corsika::units::si {
using MassDensityType = phys::units::quantity<phys::units::mass_density_d, double>;
using GrammageType = phys::units::quantity<phys::units::dimensions<-2, 1, 0>, double>;
using MomentumType = phys::units::quantity<momentum_d, double>;
- using CrossSectionType = phys::units::quantity<sigma_d, double>;
+ 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>;
} // end namespace corsika::units::si
@@ -76,11 +77,8 @@ namespace phys {
QUANTITY_DEFINE_SCALING_LITERALS(barn, corsika::units::si::sigma_d,
magnitude(corsika::units::si::constants::barn))
- QUANTITY_DEFINE_SCALING_LITERALS(meter, length_d,
- magnitude(corsika::units::si::constants::meter))
-
- QUANTITY_DEFINE_SCALING_LITERALS(newton_second, corsika::units::si::momentum_d,
- magnitude(corsika::units::si::newton_second))
+ QUANTITY_DEFINE_SCALING_LITERALS(Ns, corsika::units::si::momentum_d,
+ magnitude(1_m * 1_kg / 1_s))
} // namespace literals
} // namespace units
diff --git a/Processes/TrackingLine/TrackingLine.h b/Processes/TrackingLine/TrackingLine.h
index 205991ecf..398d4c4e5 100644
--- a/Processes/TrackingLine/TrackingLine.h
+++ b/Processes/TrackingLine/TrackingLine.h
@@ -41,14 +41,14 @@ namespace corsika::process {
public:
std::optional<std::pair<corsika::units::si::TimeType, corsika::units::si::TimeType>>
- TimeOfIntersection(corsika::geometry::Line const& traj,
+ TimeOfIntersection(corsika::geometry::Line const& line,
geometry::Sphere const& sphere) {
using namespace corsika::units::si;
auto const& cs = fEnvironment.GetCoordinateSystem();
geometry::Point const origin(cs, 0_m, 0_m, 0_m);
- auto const r0 = (traj.GetR0() - origin);
- auto const v0 = traj.GetV0();
+ auto const r0 = (line.GetR0() - origin);
+ auto const v0 = line.GetV0();
auto const c0 = (sphere.GetCenter() - origin);
auto const alpha = r0.dot(v0) - 2 * v0.dot(c0);
@@ -74,12 +74,15 @@ namespace corsika::process {
: fEnvironment(pEnv) {}
void Init() {}
- auto GetTrack(Particle& p) {
+ 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);
auto const* currentVolumeNode =
--
GitLab