diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 9c56a9a23cbdf58be5e1bd5ddd2d3736e368637c..f2326e06e981ac307784b54366d6ba001b290bb3 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -27,6 +27,8 @@
#include <corsika/units/PhysicalUnits.h>
+#include <corsika/random/RNGManager.h>
+
#include <iostream>
#include <limits>
#include <typeinfo>
@@ -116,13 +118,13 @@ public:
}
template <typename Particle>
- double MaxStepLength(Particle& p, setup::Trajectory&) const {
+ LengthType MaxStepLength(Particle& p, setup::Trajectory&) const {
const Code pid = p.GetPID();
if (isEmParticle(pid) || isInvisible(pid)) {
cout << "ProcessCut: MinStep: next cut: " << 0. << endl;
- return 0.;
+ return 0_m;
} else {
- double next_step = std::numeric_limits<double>::infinity();
+ LengthType next_step = 1_m * std::numeric_limits<double>::infinity();
cout << "ProcessCut: MinStep: next cut: " << next_step << endl;
return next_step;
}
@@ -196,8 +198,10 @@ public:
private:
};
-
int main() {
+
+ corsika::random::RNGManager::GetInstance().RegisterRandomStream("cascade");
+
corsika::environment::Environment env; // dummy environment
auto& universe = *(env.GetUniverse());
@@ -215,7 +219,8 @@ int main() {
universe.AddChild(std::move(theMedium));
- CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+ CoordinateSystem& rootCS =
+ RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
tracking_line::TrackingLine<setup::Stack> tracking(env);
stack_inspector::StackInspector<setup::Stack> p0(true);
@@ -224,9 +229,10 @@ int main() {
corsika::process::sibyll::Decay decay;
ProcessEMCut cut;
const auto sequence = /*p0 +*/ sibyll + decay + cut;
+
setup::Stack stack;
- corsika::cascade::Cascade EAS(tracking, sequence, stack);
+ corsika::cascade::Cascade EAS(env, tracking, sequence, stack);
stack.Clear();
auto particle = stack.NewParticle();
@@ -242,8 +248,9 @@ int main() {
particle.SetPosition(p);
EAS.Init();
EAS.Run();
- cout << "Result: E0=" << E0 / 1_GeV
- //<< "GeV, particles below energy threshold =" << p1.GetCount()
+ cout << "Result: E0="
+ << E0 / 1_GeV
+ //<< "GeV, particles below energy threshold =" << p1.GetCount()
<< endl;
cout << "total energy below threshold (GeV): " //<< p1.GetEnergy() / 1_GeV
<< std::endl;
diff --git a/Environment/Environment.h b/Environment/Environment.h
index cb75b64855f7d8e49341d3e8b5a17a7dc9e80b98..aa567e8764dfe981324ea53e20b0bb02b5c44a1c 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().GetRootCoordinateSystem()},
- fUniverse(std::make_unique<VolumeTreeNode<IEnvironmentModel>>(
+ : fCoordinateSystem{corsika::geometry::RootCoordinateSystem::GetInstance()
+ .GetRootCoordinateSystem()}
+ , 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 a8fa5f7283cfe0e776383f24915b98748d54349f..566b565e7ed4b7d03f37cf9a4a7c0baf7b232409 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;
+ return pTo * fDensity;
}
- corsika::units::si::TimeType FromGrammage(
+ corsika::units::si::LengthType 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 7cd9d077de0bc6f84174789771714978d8ac3f94..0edd50502ef35bdea96a608e0ab64b157394ff6d 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -12,34 +12,34 @@
#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 <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;
public:
- Cascade(Tracking& tr, ProcessList& pl, Stack& stack)
- : fTracking(tr)
- , fProcesseList(pl)
+ Cascade(corsika::environment::Environment const& env, Tracking& tr, ProcessList& pl,
+ Stack& stack)
+ : fEnvironment(env)
+ , fTracking(tr)
+ , fProcessSequence(pl)
, fStack(stack) {}
void Init() {
fTracking.Init();
- fProcesseList.Init();
+ fProcessSequence.Init();
fStack.Init();
}
@@ -56,85 +56,96 @@ 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");
+ using namespace corsika::units::si;
// determine geometric tracking
corsika::setup::Trajectory step = fTracking.GetTrack(particle);
// determine combined total interaction length (inverse)
- const double total_inv_lambda =
- fProcesseList.GetTotalInverseInteractionLength(particle, step);
+ InverseGrammageType const total_inv_lambda =
+ fProcessSequence.GetTotalInverseInteractionLength(particle, step);
+
+ // 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);
- // sample random exponential step length
- const double sample_step = rmng() / (double)rmng.max();
- const double next_interact = -log(sample_step) / total_inv_lambda;
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
- const double 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
- const double total_inv_lifetime = fProcesseList.GetTotalInverseLifetime(particle);
+ InverseTimeType const total_inv_lifetime =
+ fProcessSequence.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;
- // ....
-
// 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:
+ // 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
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
} 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 =
- fProcesseList.GetTotalInverseInteractionLength(particle, step);
- const double sample_process = rmng() / (double)rmng.max();
- double inv_lambda_count = 0;
- fProcesseList.SelectInteraction(particle, fStack, actual_inv_length,
- sample_process, inv_lambda_count);
+
+ InverseGrammageType const actual_inv_length =
+ fProcessSequence.GetTotalInverseInteractionLength(particle, step);
+
+ 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 =
- fProcesseList.GetTotalInverseLifetime(particle);
- const double sample_process = rmng() / (double)rmng.max();
- double inv_decay_count = 0;
- fProcesseList.SelectDecay(particle, fStack, actual_decay_time, sample_process,
- inv_decay_count);
+ InverseTimeType const actual_decay_time =
+ fProcessSequence.GetTotalInverseLifetime(particle);
+
+ 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);
}
}
}
@@ -142,8 +153,11 @@ namespace corsika::cascade {
private:
Tracking& fTracking;
- ProcessList& fProcesseList;
+ ProcessList& fProcessSequence;
Stack& fStack;
+ corsika::environment::Environment const& fEnvironment;
+ corsika::random::RNG& fRNG =
+ corsika::random::RNGManager::GetInstance().GetRandomStream("cascade");
};
} // namespace corsika::cascade
diff --git a/Framework/Cascade/testCascade.cc b/Framework/Cascade/testCascade.cc
index 875ff1df3a387d5665c8c7bb8b4156c1ab929690..3d7d0741948e5a9d3668501ea1c78f4398763b32 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,8 +116,9 @@ TEST_CASE("Cascade", "[Cascade]") {
const auto sequence = p0 + p1;
setup::Stack stack;
- corsika::cascade::Cascade EAS(tracking, sequence, stack);
- CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+ corsika::cascade::Cascade EAS(env, tracking, sequence, stack);
+ CoordinateSystem const& rootCS =
+ RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
stack.Clear();
auto particle = stack.NewParticle();
diff --git a/Framework/Geometry/Helix.h b/Framework/Geometry/Helix.h
index 2c8f6b8792cb204e652d1222b38f7849b70e0792..27e3575d37026e16ed401974788b80cb35b08bde 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 a79ba9ce85167ebdf5cef65622e060e844cc22b4..1b45f535a9c01acaa0084a42ca3b532360fa9b59 100644
--- a/Framework/Geometry/Line.h
+++ b/Framework/Geometry/Line.h
@@ -32,6 +32,10 @@ namespace corsika::geometry {
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 {
return v0.norm() * (t2 - t1);
diff --git a/Framework/Geometry/RootCoordinateSystem.h b/Framework/Geometry/RootCoordinateSystem.h
index a9df34edfeda18b9c7a058503ac163a0f18d4c05..3351d104d760e6d1a27269f3aa8cae1eafd1049b 100644
--- a/Framework/Geometry/RootCoordinateSystem.h
+++ b/Framework/Geometry/RootCoordinateSystem.h
@@ -22,7 +22,9 @@ namespace corsika::geometry {
public:
corsika::geometry::CoordinateSystem& GetRootCoordinateSystem() { return fRootCS; }
- const corsika::geometry::CoordinateSystem& GetRootCoordinateSystem() const { return fRootCS; }
+ const corsika::geometry::CoordinateSystem& GetRootCoordinateSystem() const {
+ return fRootCS;
+ }
private:
corsika::geometry::CoordinateSystem fRootCS; // THIS IS IT
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 bf6f8760c267856b8cd46f15590c453cc94db741..f2ac12e13c912bb36d1e5956190db9e450cb7ba9 100644
--- a/Framework/Geometry/testGeometry.cc
+++ b/Framework/Geometry/testGeometry.cc
@@ -29,7 +29,8 @@ using namespace corsika::units::si;
double constexpr absMargin = 1.0e-8;
TEST_CASE("transformations between CoordinateSystems") {
- CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+ CoordinateSystem& rootCS =
+ RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
REQUIRE(CoordinateSystem::GetTransformation(rootCS, rootCS)
.isApprox(EigenTransform::Identity()));
@@ -128,7 +129,8 @@ TEST_CASE("transformations between CoordinateSystems") {
}
TEST_CASE("Sphere") {
- CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+ CoordinateSystem& rootCS =
+ RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
Point center(rootCS, {0_m, 3_m, 4_m});
Sphere sphere(center, 5_m);
@@ -147,24 +149,34 @@ TEST_CASE("Sphere") {
}
TEST_CASE("Trajectories") {
- CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+ CoordinateSystem& rootCS =
+ RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
Point r0(rootCS, {0_m, 0_m, 0_m});
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));
+
auto const t = 1_s;
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 +186,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 +201,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/Particles/ParticleProperties.h b/Framework/Particles/ParticleProperties.h
index 91bb74eefd485b3f8afbd67fc0eaeabdbf4f7b42..db3f10fe3a411321b67c1f3fded3536a62a3816e 100644
--- a/Framework/Particles/ParticleProperties.h
+++ b/Framework/Particles/ParticleProperties.h
@@ -54,7 +54,7 @@ namespace corsika::particles {
bool constexpr IsNucleus(Code const);
int constexpr GetNucleusA(Code const);
int constexpr GetNucleusZ(Code const);
-
+
#include <corsika/particles/GeneratedParticleProperties.inc>
/*!
@@ -84,7 +84,8 @@ namespace corsika::particles {
}
corsika::units::si::TimeType constexpr GetLifetime(Code const p) {
- return detail::lifetime[static_cast<CodeIntType const>(p)] * corsika::units::si::second;
+ return detail::lifetime[static_cast<CodeIntType const>(p)] *
+ corsika::units::si::second;
}
bool constexpr IsNucleus(Code const p) {
@@ -92,14 +93,13 @@ namespace corsika::particles {
}
int constexpr GetNucleusA(Code const p) {
- return detail::nucleusA[static_cast<CodeIntType const>(p)];
+ return detail::nucleusA[static_cast<CodeIntType const>(p)];
}
int constexpr GetNucleusZ(Code const p) {
return detail::nucleusZ[static_cast<CodeIntType const>(p)];
}
-
namespace io {
std::ostream& operator<<(std::ostream& stream, Code const p);
diff --git a/Framework/Particles/testParticles.cc b/Framework/Particles/testParticles.cc
index bf307798b57d4998115a54321de965e6e53d8bb3..219e8dd8904d925eb365b76e38f6b3b9218b5f90 100644
--- a/Framework/Particles/testParticles.cc
+++ b/Framework/Particles/testParticles.cc
@@ -58,12 +58,12 @@ TEST_CASE("ParticleProperties", "[Particles]") {
REQUIRE(GetLifetime(Code::Electron) ==
std::numeric_limits<double>::infinity() * corsika::units::si::second);
REQUIRE(GetLifetime(Code::DPlus) < GetLifetime(Code::Gamma));
- REQUIRE(GetLifetime(Code::RhoPlus)/corsika::units::si::second ==
- (Approx(4.414566727909413e-24).epsilon(1e-3)));
- REQUIRE(GetLifetime(Code::SigmaMinusBar)/corsika::units::si::second ==
- (Approx(8.018880848563575e-11).epsilon(1e-5)));
- REQUIRE(GetLifetime(Code::MuPlus)/corsika::units::si::second ==
- (Approx(2.1970332555864364e-06).epsilon(1e-5)));
+ REQUIRE(GetLifetime(Code::RhoPlus) / corsika::units::si::second ==
+ (Approx(4.414566727909413e-24).epsilon(1e-3)));
+ REQUIRE(GetLifetime(Code::SigmaMinusBar) / corsika::units::si::second ==
+ (Approx(8.018880848563575e-11).epsilon(1e-5)));
+ REQUIRE(GetLifetime(Code::MuPlus) / corsika::units::si::second ==
+ (Approx(2.1970332555864364e-06).epsilon(1e-5)));
}
SECTION("Nuclei") {
@@ -78,8 +78,5 @@ TEST_CASE("ParticleProperties", "[Particles]") {
REQUIRE(GetNucleusA(Code::Tritium) == 3);
REQUIRE(Hydrogen::GetNucleusZ() == 1);
REQUIRE(Tritium::GetNucleusA() == 3);
-
-
}
-
}
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 014630e1a46c076f5f8347686f1d047771c0b29f..d58c14b5083cfda5db73f468357fa28e38122248 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>
@@ -160,7 +161,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,32 +175,41 @@ namespace corsika::process {
}
template <typename Particle, typename Track>
- inline double MaxStepLength(Particle& p, Track& track) const {
- double max_length = std::numeric_limits<double>::infinity();
+ 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>
- inline double 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>
- inline double 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>
- inline double GetInverseInteractionLength(Particle& p, Track& t) const {
- double 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,21 @@ 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 {
+ 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_count / lambda_inv_tot) {
+ if (lambda_select < lambda_inv_count) {
A.DoInteraction(p, s);
return EProcessReturn::eInteracted;
}
@@ -235,14 +245,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 +261,21 @@ namespace corsika::process {
}
template <typename Particle>
- inline double GetTotalLifetime(Particle& p) const {
+ corsika::units::si::TimeType GetTotalLifetime(Particle& p) const {
return 1. / GetInverseLifetime(p);
}
template <typename Particle>
- inline double GetTotalInverseLifetime(Particle& p) const {
+ corsika::units::si::InverseTimeType GetTotalInverseLifetime(Particle& p) const {
return GetInverseLifetime(p);
}
template <typename Particle>
- inline double GetInverseLifetime(Particle& p) const {
- double 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 +289,20 @@ 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, 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_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 +310,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;
}
@@ -314,7 +326,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();
}
@@ -349,7 +361,7 @@ namespace corsika::process {
OPSEQ(DecayProcess, DecayProcess)
template <typename A, typename B>
- struct is_process_sequence<corsika::process::ProcessSequence<A, B> > {
+ struct is_process_sequence<corsika::process::ProcessSequence<A, B> > {
static const bool value = true;
};
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/Random/CMakeLists.txt b/Framework/Random/CMakeLists.txt
index 3c1947d2508b318eafe6ab426bceac9af7b9e752..230f84f67f4ed6eb02652c5cb10853955f7de776 100644
--- a/Framework/Random/CMakeLists.txt
+++ b/Framework/Random/CMakeLists.txt
@@ -7,6 +7,7 @@ set (
set (
CORSIKArandom_HEADERS
RNGManager.h
+ UniformRealDistribution.h
)
set (
diff --git a/Framework/Random/UniformRealDistribution.h b/Framework/Random/UniformRealDistribution.h
new file mode 100644
index 0000000000000000000000000000000000000000..c90ab674ea41976647320a64095a5e8b68d6484e
--- /dev/null
+++ b/Framework/Random/UniformRealDistribution.h
@@ -0,0 +1,32 @@
+#ifndef _include_random_distributions_h
+#define _include_random_distributions_h
+
+#include <corsika/units/PhysicalUnits.h>
+#include <random>
+
+namespace corsika::random {
+
+ template <class TQuantity>
+ class UniformRealDistribution {
+ using RealType = typename TQuantity::value_type;
+ std::uniform_real_distribution<RealType> dist{RealType(0.), RealType(1.)};
+
+ TQuantity const a, b;
+
+ public:
+ UniformRealDistribution(TQuantity b)
+ : a{TQuantity(phys::units::detail::magnitude_tag, 0)}
+ , b(b) {}
+ UniformRealDistribution(TQuantity a, TQuantity b)
+ : a(a)
+ , b(b) {}
+
+ template <class Generator>
+ TQuantity operator()(Generator& g) {
+ return a + dist(g) * (b - a);
+ }
+ };
+
+} // namespace corsika::random
+
+#endif
diff --git a/Framework/Random/testRandom.cc b/Framework/Random/testRandom.cc
index 041a96c21ff241e3a86a903a887e521ff17c9e29..56d8157a74530d404f5d8ba1c8cf8aefe9f75f99 100644
--- a/Framework/Random/testRandom.cc
+++ b/Framework/Random/testRandom.cc
@@ -1,4 +1,3 @@
-
/**
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
@@ -14,7 +13,11 @@
#include <catch2/catch.hpp>
#include <corsika/random/RNGManager.h>
+#include <corsika/random/UniformRealDistribution.h>
+#include <corsika/units/PhysicalUnits.h>
#include <iostream>
+#include <limits>
+#include <random>
using namespace corsika::random;
@@ -37,3 +40,46 @@ SCENARIO("random-number streams can be registered and retrieved") {
}
}
}
+
+TEST_CASE("UniformRealDistribution") {
+ using namespace corsika::units::si;
+ std::mt19937 rng;
+
+ corsika::random::UniformRealDistribution<LengthType> dist(1_m, 2_m);
+
+ SECTION("range") {
+ corsika::random::UniformRealDistribution<LengthType> dist(1_m, 2_m);
+
+ LengthType min =
+ +1_m * std::numeric_limits<typename LengthType::value_type>::infinity();
+ LengthType max =
+ -1_m * std::numeric_limits<typename LengthType::value_type>::infinity();
+
+ for (int i{0}; i < 1'000'000; ++i) {
+ LengthType x = dist(rng);
+ min = std::min(min, x);
+ max = std::max(max, x);
+ }
+
+ CHECK(min / 1_m == Approx(1.));
+ CHECK(max / 2_m == Approx(1.));
+ }
+
+ SECTION("range") {
+ corsika::random::UniformRealDistribution<LengthType> dist(18_cm);
+
+ LengthType min =
+ +1_m * std::numeric_limits<typename LengthType::value_type>::infinity();
+ LengthType max =
+ -1_m * std::numeric_limits<typename LengthType::value_type>::infinity();
+
+ for (int i{0}; i < 1'000'000; ++i) {
+ LengthType x = dist(rng);
+ min = std::min(min, x);
+ max = std::max(max, x);
+ }
+
+ CHECK(min / 1_m == Approx(0.).margin(1e-3));
+ CHECK(max / 18_cm == Approx(1.));
+ }
+}
diff --git a/Framework/Units/PhysicalUnits.h b/Framework/Units/PhysicalUnits.h
index 7083ce743811176e2c807f2d421dde2273fbb673..03029541311ccbbe259a2f49901e96ac76a4cdb9 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,8 +53,13 @@ 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>;
+ using InverseGrammageType =
+ phys::units::quantity<phys::units::dimensions<2, -1, 0>, double>;
} // end namespace corsika::units::si
/**
@@ -73,14 +77,11 @@ 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(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/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/NullModel/NullModel.cc b/Processes/NullModel/NullModel.cc
index 631891e63d11cca0fc7485c4e549496037e0b19e..313d38034d3c90925eb461aa59b1399af8308e4d 100644
--- a/Processes/NullModel/NullModel.cc
+++ b/Processes/NullModel/NullModel.cc
@@ -11,7 +11,6 @@
#include <corsika/process/null_model/NullModel.h>
-
#include <corsika/logging/Logger.h>
#include <corsika/setup/SetupTrajectory.h>
@@ -31,7 +30,7 @@ NullModel<Stack>::~NullModel() {}
template <typename Stack>
process::EProcessReturn NullModel<Stack>::DoContinuous(Particle&, setup::Trajectory&,
- Stack& ) const {
+ Stack&) const {
return EProcessReturn::eOk;
}
@@ -41,10 +40,8 @@ double NullModel<Stack>::MaxStepLength(Particle&, setup::Trajectory&) const {
}
template <typename Stack>
-void NullModel<Stack>::Init() {
-}
+void NullModel<Stack>::Init() {}
#include <corsika/setup/SetupStack.h>
template class process::null_model::NullModel<setup::Stack>;
-
diff --git a/Processes/NullModel/NullModel.h b/Processes/NullModel/NullModel.h
index e7eac1bb0196c77f0ced559ac6145455615fd00b..367c907f38df5ab6b04bbc23ccaefc24a54df4ad 100644
--- a/Processes/NullModel/NullModel.h
+++ b/Processes/NullModel/NullModel.h
@@ -31,7 +31,6 @@ namespace corsika::process {
void Init();
EProcessReturn DoContinuous(Particle&, corsika::setup::Trajectory&, Stack& s) const;
double MaxStepLength(Particle&, corsika::setup::Trajectory&) const;
-
};
} // namespace null_model
diff --git a/Processes/NullModel/testNullModel.cc b/Processes/NullModel/testNullModel.cc
index 908d7370224a4237bb1332e6b6149176eace39ee..ffcba5e0351acefa519f6a28e003be7999741f68 100644
--- a/Processes/NullModel/testNullModel.cc
+++ b/Processes/NullModel/testNullModel.cc
@@ -16,8 +16,8 @@
#include <corsika/process/null_model/NullModel.h>
#include <corsika/geometry/Point.h>
-#include <corsika/geometry/Vector.h>
#include <corsika/geometry/RootCoordinateSystem.h>
+#include <corsika/geometry/Vector.h>
#include <corsika/units/PhysicalUnits.h>
@@ -30,24 +30,24 @@ using namespace corsika;
TEST_CASE("NullModel", "[processes]") {
- auto const& cs = geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
- geometry::Point const origin(cs, {0_m, 0_m, 0_m});
- geometry::Vector<corsika::units::si::SpeedType::dimension_type> v(cs, 0_m / second,
- 0_m / second, 1_m / second);
+ auto const& cs =
+ geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+ geometry::Point const origin(cs, {0_m, 0_m, 0_m});
+ geometry::Vector<corsika::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);
setup::Stack stack;
auto particle = stack.NewParticle();
-
+
SECTION("interface") {
NullModel<setup::Stack> model;
model.Init();
- [[maybe_unused]] const process::EProcessReturn ret = model.DoContinuous(particle, track, stack);
+ [[maybe_unused]] const process::EProcessReturn ret =
+ model.DoContinuous(particle, track, stack);
[[maybe_unused]] const double length = model.MaxStepLength(particle, track);
-
-
}
}
diff --git a/Processes/Sibyll/Decay.h b/Processes/Sibyll/Decay.h
index fd1a8f4237c63d7531e07bd865da4de6e1191b53..050466793b518aed2e07bed12a356b9933ba30cf 100644
--- a/Processes/Sibyll/Decay.h
+++ b/Processes/Sibyll/Decay.h
@@ -1,12 +1,11 @@
#ifndef _include_corsika_process_sibyll_decay_h_
#define _include_corsika_process_sibyll_decay_h_
-#include <corsika/process/sibyll/SibStack.h>
-#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/process/DecayProcess.h>
-
-#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/process/sibyll/ParticleConversion.h>
+#include <corsika/process/sibyll/SibStack.h>
#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
#include <corsika/particles/ParticleProperties.h>
@@ -59,8 +58,8 @@ namespace corsika::process {
void setTrackedParticlesStable() {
/*
- Sibyll is hadronic generator
- only hadrons decay
+ Sibyll is hadronic generator
+ only hadrons decay
*/
// set particles unstable
setHadronsUnstable();
@@ -73,12 +72,12 @@ namespace corsika::process {
ds.NewParticle().SetPID(particles::Code::KMinus);
ds.NewParticle().SetPID(particles::Code::K0Long);
ds.NewParticle().SetPID(particles::Code::K0Short);
-
+
for (auto& p : ds) {
- int s_id = process::sibyll::ConvertToSibyllRaw(p.GetPID());
- // set particle stable by setting table value negative
- s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]);
- p.Delete();
+ int s_id = process::sibyll::ConvertToSibyllRaw(p.GetPID());
+ // set particle stable by setting table value negative
+ s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]);
+ p.Delete();
}
}
@@ -86,7 +85,7 @@ namespace corsika::process {
public:
Decay() {}
void Init() {
- setHadronsUnstable();
+ setHadronsUnstable();
setTrackedParticlesStable();
}
@@ -112,11 +111,11 @@ namespace corsika::process {
friend void setHadronsUnstable();
template <typename Particle>
- double GetLifetime(Particle& p) const {
+ corsika::units::si::TimeType GetLifetime(Particle& p) const {
corsika::units::hep::EnergyType E = p.GetEnergy();
corsika::units::hep::MassType m = corsika::particles::GetMass(p.GetPID());
- //const MassDensityType density = 1.25e-3 * kilogram / (1_cm * 1_cm * 1_cm);
+ // const MassDensityType density = 1.25e-3 * kilogram / (1_cm * 1_cm * 1_cm);
const double gamma = E / m;
@@ -128,11 +127,11 @@ namespace corsika::process {
// return as column density
// const double x0 = density * t0 * gamma * constants::c / kilogram * 1_cm * 1_cm;
// cout << "Decay: MinStep: x0: " << x0 << endl;
- const double lifetime = gamma*t0 / 1_s;
- cout << "Decay: MinStep: tau: " << lifetime << endl;
- //int a = 1;
- //const double x = -x0 * log(s_rndm_(a));
- //cout << "Decay: next decay: " << x << endl;
+ corsika::units::si::TimeType const lifetime = gamma * t0;
+ cout << "Decay: MinStep: tau: " << lifetime << endl;
+ // int a = 1;
+ // const double x = -x0 * log(s_rndm_(a));
+ // cout << "Decay: next decay: " << x << endl;
return lifetime;
}
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
index fa22f7fc71c18705558cab61987dd722fb8aba2d..24401ad7c66867370c18ce4101100d998d62b24e 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/Processes/Sibyll/Interaction.h
@@ -6,272 +6,256 @@
//#include <corsika/setup/SetupStack.h>
//#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/process/sibyll/SibStack.h>
#include <corsika/process/sibyll/sibyll2.3c.h>
-#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/particles/ParticleProperties.h>
-#include <corsika/units/PhysicalUnits.h>
#include <corsika/random/RNGManager.h>
+#include <corsika/units/PhysicalUnits.h>
using namespace corsika;
using namespace corsika::process::sibyll;
using namespace corsika::units::si;
-
namespace corsika::process::sibyll {
- class Interaction : public corsika::process::InteractionProcess<Interaction> {
+ class Interaction : public corsika::process::InteractionProcess<Interaction> {
public:
-
Interaction() {}
~Interaction() {}
-
+
void Init() {
corsika::random::RNGManager& rmng = corsika::random::RNGManager::GetInstance();
rmng.RegisterRandomStream("s_rndm");
-
+
// test random number generator
std::cout << "Interaction: "
- << " test sequence of random numbers." << std::endl;
+ << " test sequence of random numbers." << std::endl;
int a = 0;
for (int i = 0; i < 8; ++i) std::cout << i << " " << s_rndm_(a) << std::endl;
-
+
// initialize Sibyll
sibyll_ini_();
}
- //void setTrackedParticlesStable();
+ // void setTrackedParticlesStable();
template <typename Particle, typename Track>
- double GetInteractionLength(Particle& p, Track&) const {
-
- // coordinate system, get global frame of reference
- CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
-
- const particles::Code corsikaBeamId = p.GetPID();
-
- // beam particles for sibyll : 1, 2, 3 for p, pi, k
- // read from cross section code table
- int kBeam = process::sibyll::GetSibyllXSCode(corsikaBeamId);
-
- bool kInteraction = process::sibyll::CanInteract(corsikaBeamId);
-
- /*
- the target should be defined by the Environment,
- ideally as full particle object so that the four momenta
- and the boosts can be defined..
- */
- // target nuclei: A < 18
- // FOR NOW: assume target is oxygen
- int kTarget = 16;
-
- EnergyType Etot = p.GetEnergy() + kTarget * corsika::particles::Proton::GetMass();
- super_stupid::MomentumVector Ptot(
- rootCS, {0.0_newton_second, 0.0_newton_second, 0.0_newton_second});
- // FOR NOW: assume target is at rest
- super_stupid::MomentumVector pTarget(
- rootCS, {0.0_newton_second, 0.0_newton_second, 0.0_newton_second});
- Ptot += p.GetMomentum();
- Ptot += pTarget;
- // calculate cm. energy
- EnergyType sqs = sqrt(Etot * Etot - Ptot.squaredNorm() * constants::cSquared);
- double Ecm = sqs / 1_GeV;
-
- std::cout << "Interaction: "
- << "MinStep: input en: " << p.GetEnergy() / 1_GeV << endl
- << " beam can interact:" << kBeam << endl
- << " beam XS code:" << kBeam << endl
- << " beam pid:" << p.GetPID() << endl
- << " target mass number:" << kTarget << std::endl;
-
- double int_length = 0;
- if (kInteraction) {
-
- double prodCrossSection, dummy, dum1, dum2, dum3, dum4;
- double dumdif[3];
-
- if (kTarget == 1)
- sib_sigma_hp_(kBeam, Ecm, dum1, dum2, prodCrossSection, dumdif, dum3, dum4);
- else
- sib_sigma_hnuc_(kBeam, kTarget, Ecm, prodCrossSection, dummy);
-
- std::cout << "Interaction: "
- << "MinStep: sibyll return: " << prodCrossSection << std::endl;
- CrossSectionType sig = prodCrossSection * 1_mbarn;
- std::cout << "Interaction: "
- << "MinStep: CrossSection (mb): " << sig / 1_mbarn << std::endl;
+ corsika::units::si::GrammageType GetInteractionLength(Particle& p, Track&) const {
+ // coordinate system, get global frame of reference
+ CoordinateSystem& rootCS =
+ RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
- const MassType nucleon_mass = 0.93827_GeV / corsika::units::si::constants::cSquared;
- std::cout << "Interaction: "
- << "nucleon mass " << nucleon_mass << std::endl;
- // calculate interaction length in medium
- int_length = kTarget * (nucleon_mass / 1_g) / (sig / 1_cmeter / 1_cmeter);
- // pick random step lenth
- std::cout << "Interaction: "
- << "interaction length (g/cm2): " << int_length << std::endl;
- } else
- int_length = std::numeric_limits<double>::infinity();
-
- /*
- what are the units of the output? slant depth or 3space length?
-
- */
- return int_length;
- //
- // int a = 0;
- // const double next_step = -int_length * log(s_rndm_(a));
- // std::cout << "Interaction: "
- // << "next step (g/cm2): " << next_step << std::endl;
- // return next_step;
- }
-
+ const particles::Code corsikaBeamId = p.GetPID();
- template <typename Particle, typename Stack>
- corsika::process::EProcessReturn DoInteraction(Particle& p, Stack& s) const {
- cout << "ProcessSibyll: "
- << "DoInteraction: " << p.GetPID() << " interaction? "
- << process::sibyll::CanInteract(p.GetPID()) << endl;
- if (process::sibyll::CanInteract(p.GetPID())) {
- cout << "defining coordinates" << endl;
- // coordinate system, get global frame of reference
- CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+ // beam particles for sibyll : 1, 2, 3 for p, pi, k
+ // read from cross section code table
+ int kBeam = process::sibyll::GetSibyllXSCode(corsikaBeamId);
- QuantityVector<length_d> const coordinates{0_m, 0_m, 0_m};
- Point pOrig(rootCS, coordinates);
+ bool kInteraction = process::sibyll::CanInteract(corsikaBeamId);
/*
the target should be defined by the Environment,
ideally as full particle object so that the four momenta
and the boosts can be defined..
+ */
+ // target nuclei: A < 18
+ // FOR NOW: assume target is oxygen
+ int kTarget = 16;
+
+ EnergyType Etot = p.GetEnergy() + kTarget * corsika::particles::Proton::GetMass();
+ super_stupid::MomentumVector Ptot(rootCS, {0.0_Ns, 0.0_Ns, 0.0_Ns});
+ // FOR NOW: assume target is at rest
+ super_stupid::MomentumVector pTarget(rootCS, {0.0_Ns, 0.0_Ns, 0.0_Ns});
+ Ptot += p.GetMomentum();
+ Ptot += pTarget;
+ // calculate cm. energy
+ EnergyType sqs = sqrt(Etot * Etot - Ptot.squaredNorm() * constants::cSquared);
+ double Ecm = sqs / 1_GeV;
- here we need: GetTargetMassNumber() or GetTargetPID()??
- GetTargetMomentum() (zero in EAS)
- */
- // FOR NOW: set target to proton
- int kTarget = 1; // env.GetTargetParticle().GetPID();
-
- cout << "defining target momentum.." << endl;
- // FOR NOW: target is always at rest
- const EnergyType Etarget = 0. * 1_GeV + corsika::particles::Proton::GetMass();
- const auto pTarget = super_stupid::MomentumVector(
- rootCS, 0. * 1_GeV / constants::c, 0. * 1_GeV / constants::c,
- 0. * 1_GeV / constants::c);
- cout << "target momentum (GeV/c): "
- << pTarget.GetComponents() / 1_GeV * constants::c << endl;
- cout << "beam momentum (GeV/c): "
- << p.GetMomentum().GetComponents() / 1_GeV * constants::c << endl;
-
- // get energy of particle from stack
- /*
- stack is in GeV in lab. frame
- convert to GeV in cm. frame
- (assuming proton at rest as target AND
- assuming no pT, i.e. shower frame-z is aligned with hadron-int-frame-z)
- */
- // total energy: E_beam + E_target
- // in lab. frame: E_beam + m_target*c**2
- EnergyType E = p.GetEnergy();
- EnergyType Etot = E + Etarget;
- // total momentum
- super_stupid::MomentumVector Ptot = p.GetMomentum(); // + pTarget;
- // invariant mass, i.e. cm. energy
- EnergyType Ecm = sqrt(Etot * Etot -
- Ptot.squaredNorm() *
- constants::cSquared); // sqrt( 2. * E * 0.93827_GeV );
- /*
- get transformation between Stack-frame and SibStack-frame
- for EAS Stack-frame is lab. frame, could be different for CRMC-mode
- the transformation should be derived from the input momenta
- */
- const double gamma = Etot / Ecm;
- const auto gambet = Ptot / (Ecm / constants::c);
-
- std::cout << "Interaction: "
- << " DoDiscrete: gamma:" << gamma << endl;
std::cout << "Interaction: "
- << " DoDiscrete: gambet:" << gambet.GetComponents() << endl;
+ << "MinStep: input en: " << p.GetEnergy() / 1_GeV << endl
+ << " beam can interact:" << kBeam << endl
+ << " beam XS code:" << kBeam << endl
+ << " beam pid:" << p.GetPID() << endl
+ << " target mass number:" << kTarget << std::endl;
- int kBeam = process::sibyll::ConvertToSibyllRaw(p.GetPID());
+ if (kInteraction) {
- std::cout << "Interaction: "
- << " DoInteraction: E(GeV):" << E / 1_GeV << " Ecm(GeV): " << Ecm / 1_GeV
- << std::endl;
- if (E < 8.5_GeV || Ecm < 10_GeV) {
+ double prodCrossSection, dummy, dum1, dum2, dum3, dum4;
+ double dumdif[3];
+
+ if (kTarget == 1)
+ sib_sigma_hp_(kBeam, Ecm, dum1, dum2, prodCrossSection, dumdif, dum3, dum4);
+ else
+ sib_sigma_hnuc_(kBeam, kTarget, Ecm, prodCrossSection, dummy);
+
+ std::cout << "Interaction: "
+ << "MinStep: sibyll return: " << prodCrossSection << std::endl;
+ CrossSectionType sig = prodCrossSection * 1_mbarn;
+ std::cout << "Interaction: "
+ << "MinStep: CrossSection (mb): " << sig / 1_mbarn << std::endl;
+
+ const MassType nucleon_mass =
+ 0.93827_GeV / corsika::units::si::constants::cSquared;
std::cout << "Interaction: "
- << " DoInteraction: dropping particle.." << std::endl;
- p.Delete();
+ << "nucleon mass " << nucleon_mass << std::endl;
+ // calculate interaction length in medium
+ GrammageType int_length = kTarget * nucleon_mass / sig;
+ // pick random step lenth
+ std::cout << "Interaction: "
+ << "interaction length (g/cm2): " << int_length << std::endl;
+
+ return int_length;
} else {
- // Sibyll does not know about units..
- double sqs = Ecm / 1_GeV;
- // running sibyll, filling stack
- sibyll_(kBeam, kTarget, sqs);
- // running decays
- //setTrackedParticlesStable();
- decsib_();
- // print final state
- int print_unit = 6;
- sib_list_(print_unit);
-
- // delete current particle
- p.Delete();
-
- // add particles from sibyll to stack
- // link to sibyll stack
- SibStack ss;
-
- // SibStack does not know about momentum yet so we need counter to access momentum
- // array in Sibyll
- int i = -1;
- super_stupid::MomentumVector Ptot_final(
- rootCS, {0.0_newton_second, 0.0_newton_second, 0.0_newton_second});
- for (auto& psib : ss) {
- ++i;
- // skip particles that have decayed in Sibyll
- if (abs(s_plist_.llist[i]) > 100) continue;
-
- // transform energy to lab. frame, primitve
- // compute beta_vec * p_vec
- // arbitrary Lorentz transformation based on sibyll routines
- const auto gammaBetaComponents = gambet.GetComponents();
- const auto pSibyllComponents = psib.GetMomentum().GetComponents();
- EnergyType en_lab = 0. * 1_GeV;
- MomentumType p_lab_components[3];
- en_lab = psib.GetEnergy() * gamma;
- EnergyType pnorm = 0. * 1_GeV;
- for (int j = 0; j < 3; ++j)
- pnorm += (pSibyllComponents[j] * gammaBetaComponents[j] * constants::c) /
- (gamma + 1.);
- pnorm += psib.GetEnergy();
-
- for (int j = 0; j < 3; ++j) {
- p_lab_components[j] = pSibyllComponents[j] - (-1) * pnorm *
- gammaBetaComponents[j] /
- constants::c;
- en_lab -=
- (-1) * pSibyllComponents[j] * gammaBetaComponents[j] * constants::c;
- }
+ return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
+ }
+ }
- // add to corsika stack
- auto pnew = s.NewParticle();
- pnew.SetEnergy(en_lab);
- pnew.SetPID(process::sibyll::ConvertFromSibyll(psib.GetPID()));
+ template <typename Particle, typename Stack>
+ corsika::process::EProcessReturn DoInteraction(Particle& p, Stack& s) const {
+ cout << "ProcessSibyll: "
+ << "DoInteraction: " << p.GetPID() << " interaction? "
+ << process::sibyll::CanInteract(p.GetPID()) << endl;
+ if (process::sibyll::CanInteract(p.GetPID())) {
+ cout << "defining coordinates" << endl;
+ // coordinate system, get global frame of reference
+ CoordinateSystem& rootCS =
+ RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+
+ QuantityVector<length_d> const coordinates{0_m, 0_m, 0_m};
+ Point pOrig(rootCS, coordinates);
+
+ /*
+ the target should be defined by the Environment,
+ ideally as full particle object so that the four momenta
+ and the boosts can be defined..
+
+ here we need: GetTargetMassNumber() or GetTargetPID()??
+ GetTargetMomentum() (zero in EAS)
+ */
+ // FOR NOW: set target to proton
+ int kTarget = 1; // env.GetTargetParticle().GetPID();
+
+ cout << "defining target momentum.." << endl;
+ // FOR NOW: target is always at rest
+ const EnergyType Etarget = 0. * 1_GeV + corsika::particles::Proton::GetMass();
+ const auto pTarget = super_stupid::MomentumVector(
+ rootCS, 0. * 1_GeV / constants::c, 0. * 1_GeV / constants::c,
+ 0. * 1_GeV / constants::c);
+ cout << "target momentum (GeV/c): "
+ << pTarget.GetComponents() / 1_GeV * constants::c << endl;
+ cout << "beam momentum (GeV/c): "
+ << p.GetMomentum().GetComponents() / 1_GeV * constants::c << endl;
+
+ // get energy of particle from stack
+ /*
+ stack is in GeV in lab. frame
+ convert to GeV in cm. frame
+ (assuming proton at rest as target AND
+ assuming no pT, i.e. shower frame-z is aligned with hadron-int-frame-z)
+ */
+ // total energy: E_beam + E_target
+ // in lab. frame: E_beam + m_target*c**2
+ EnergyType E = p.GetEnergy();
+ EnergyType Etot = E + Etarget;
+ // total momentum
+ super_stupid::MomentumVector Ptot = p.GetMomentum(); // + pTarget;
+ // invariant mass, i.e. cm. energy
+ EnergyType Ecm =
+ sqrt(Etot * Etot - Ptot.squaredNorm() *
+ constants::cSquared); // sqrt( 2. * E * 0.93827_GeV );
+ /*
+ get transformation between Stack-frame and SibStack-frame
+ for EAS Stack-frame is lab. frame, could be different for CRMC-mode
+ the transformation should be derived from the input momenta
+ */
+ const double gamma = Etot / Ecm;
+ const auto gambet = Ptot / (Ecm / constants::c);
+
+ std::cout << "Interaction: "
+ << " DoDiscrete: gamma:" << gamma << endl;
+ std::cout << "Interaction: "
+ << " DoDiscrete: gambet:" << gambet.GetComponents() << endl;
+
+ int kBeam = process::sibyll::ConvertToSibyllRaw(p.GetPID());
- corsika::geometry::QuantityVector<momentum_d> p_lab_c{
- p_lab_components[0], p_lab_components[1], p_lab_components[2]};
- pnew.SetMomentum(super_stupid::MomentumVector(rootCS, p_lab_c));
- Ptot_final += pnew.GetMomentum();
+ std::cout << "Interaction: "
+ << " DoInteraction: E(GeV):" << E / 1_GeV
+ << " Ecm(GeV): " << Ecm / 1_GeV << std::endl;
+ if (E < 8.5_GeV || Ecm < 10_GeV) {
+ std::cout << "Interaction: "
+ << " DoInteraction: dropping particle.." << std::endl;
+ p.Delete();
+ } else {
+ // Sibyll does not know about units..
+ double sqs = Ecm / 1_GeV;
+ // running sibyll, filling stack
+ sibyll_(kBeam, kTarget, sqs);
+ // running decays
+ // setTrackedParticlesStable();
+ decsib_();
+ // print final state
+ int print_unit = 6;
+ sib_list_(print_unit);
+
+ // delete current particle
+ p.Delete();
+
+ // add particles from sibyll to stack
+ // link to sibyll stack
+ SibStack ss;
+
+ // SibStack does not know about momentum yet so we need counter to access
+ // momentum array in Sibyll
+ int i = -1;
+ super_stupid::MomentumVector Ptot_final(rootCS, {0.0_Ns, 0.0_Ns, 0.0_Ns});
+ for (auto& psib : ss) {
+ ++i;
+ // skip particles that have decayed in Sibyll
+ if (abs(s_plist_.llist[i]) > 100) continue;
+
+ // transform energy to lab. frame, primitve
+ // compute beta_vec * p_vec
+ // arbitrary Lorentz transformation based on sibyll routines
+ const auto gammaBetaComponents = gambet.GetComponents();
+ const auto pSibyllComponents = psib.GetMomentum().GetComponents();
+ EnergyType en_lab = 0. * 1_GeV;
+ MomentumType p_lab_components[3];
+ en_lab = psib.GetEnergy() * gamma;
+ EnergyType pnorm = 0. * 1_GeV;
+ for (int j = 0; j < 3; ++j)
+ pnorm += (pSibyllComponents[j] * gammaBetaComponents[j] * constants::c) /
+ (gamma + 1.);
+ pnorm += psib.GetEnergy();
+
+ for (int j = 0; j < 3; ++j) {
+ p_lab_components[j] = pSibyllComponents[j] -
+ (-1) * pnorm * gammaBetaComponents[j] / constants::c;
+ en_lab -=
+ (-1) * pSibyllComponents[j] * gammaBetaComponents[j] * constants::c;
+ }
+
+ // add to corsika stack
+ auto pnew = s.NewParticle();
+ pnew.SetEnergy(en_lab);
+ pnew.SetPID(process::sibyll::ConvertFromSibyll(psib.GetPID()));
+
+ corsika::geometry::QuantityVector<momentum_d> p_lab_c{
+ p_lab_components[0], p_lab_components[1], p_lab_components[2]};
+ pnew.SetMomentum(super_stupid::MomentumVector(rootCS, p_lab_c));
+ Ptot_final += pnew.GetMomentum();
+ }
+ // cout << "tot. momentum final (GeV/c): " << Ptot_final.GetComponents() / 1_GeV
+ // * constants::c << endl;
}
- // cout << "tot. momentum final (GeV/c): " << Ptot_final.GetComponents() / 1_GeV *
- // constants::c << endl;
}
+ return process::EProcessReturn::eOk;
}
- return process::EProcessReturn::eOk;
- }
-
};
-
-}
+
+} // namespace corsika::process::sibyll
#endif
diff --git a/Processes/Sibyll/ParticleConversion.h b/Processes/Sibyll/ParticleConversion.h
index e5e76c3b398e2d93a7af8b8fbe652dacf37ea8cb..cd96f83d5af25e4be078d693458ad3239d741aaf 100644
--- a/Processes/Sibyll/ParticleConversion.h
+++ b/Processes/Sibyll/ParticleConversion.h
@@ -25,7 +25,7 @@ namespace corsika::process::sibyll {
#include <corsika/process/sibyll/Generated.inc>
- bool constexpr KnownBySibyll(corsika::particles::Code pCode) {
+ bool constexpr KnownBySibyll(corsika::particles::Code pCode) {
return isKnown[static_cast<corsika::particles::CodeIntType>(pCode)];
}
diff --git a/Processes/Sibyll/SibStack.h b/Processes/Sibyll/SibStack.h
index 09e97b391f82a5631c7e2111f249b60fe34ef7ed..264cff1b8b97d6b3335bd22f7131b53984ae2ba6 100644
--- a/Processes/Sibyll/SibStack.h
+++ b/Processes/Sibyll/SibStack.h
@@ -1,8 +1,8 @@
#ifndef _include_sibstack_h_
#define _include_sibstack_h_
-#include <corsika/process/sibyll/sibyll2.3c.h>
#include <corsika/process/sibyll/ParticleConversion.h>
+#include <corsika/process/sibyll/sibyll2.3c.h>
#include <corsika/stack/Stack.h>
#include <corsika/units/PhysicalUnits.h>
@@ -14,75 +14,76 @@ using namespace corsika::units::si;
using namespace corsika::geometry;
namespace corsika::process::sibyll {
-
-class SibStackData {
-public:
- void Init();
+ class SibStackData {
+
+ public:
+ void Init();
- void Clear() { s_plist_.np = 0; }
+ void Clear() { s_plist_.np = 0; }
- int GetSize() const { return s_plist_.np; }
+ int GetSize() const { return s_plist_.np; }
#warning check actual capacity of sibyll stack
- int GetCapacity() const { return 8000; }
-
- void SetId(const int i, const int v) { s_plist_.llist[i] = v; }
- void SetEnergy(const int i, const EnergyType v) { s_plist_.p[3][i] = v / 1_GeV; }
- void SetMomentum(const int i, const super_stupid::MomentumVector& v) {
- auto tmp = v.GetComponents();
- for (int idx = 0; idx < 3; ++idx)
- s_plist_.p[idx][i] = tmp[idx] / 1_GeV * constants::c;
- }
-
- int GetId(const int i) const { return s_plist_.llist[i]; }
-
- EnergyType GetEnergy(const int i) const { return s_plist_.p[3][i] * 1_GeV; }
-
- super_stupid::MomentumVector GetMomentum(const int i) const {
- CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
- corsika::geometry::QuantityVector<momentum_d> components{
- s_plist_.p[0][i] * 1_GeV / constants::c,
- s_plist_.p[1][i] * 1_GeV / constants::c,
- s_plist_.p[2][i] * 1_GeV / constants::c};
- super_stupid::MomentumVector v1(rootCS, components);
- return v1;
- }
-
- void Copy(const int i1, const int i2) {
- s_plist_.llist[i1] = s_plist_.llist[i2];
- s_plist_.p[3][i1] = s_plist_.p[3][i2];
- }
-
-protected:
- void IncrementSize() { s_plist_.np++; }
- void DecrementSize() {
- if (s_plist_.np > 0) { s_plist_.np--; }
- }
-};
-
-template <typename StackIteratorInterface>
-class ParticleInterface : public ParticleBase<StackIteratorInterface> {
- using ParticleBase<StackIteratorInterface>::GetStackData;
- using ParticleBase<StackIteratorInterface>::GetIndex;
-
-public:
- void SetEnergy(const EnergyType v) { GetStackData().SetEnergy(GetIndex(), v); }
- EnergyType GetEnergy() const { return GetStackData().GetEnergy(GetIndex()); }
- void SetPID(const int v) { GetStackData().SetId(GetIndex(), v); }
- corsika::process::sibyll::SibyllCode GetPID() const {
- return static_cast<corsika::process::sibyll::SibyllCode>(
- GetStackData().GetId(GetIndex()));
- }
- super_stupid::MomentumVector GetMomentum() const {
- return GetStackData().GetMomentum(GetIndex());
- }
- void SetMomentum(const super_stupid::MomentumVector& v) {
- GetStackData().SetMomentum(GetIndex(), v);
- }
-};
-
-typedef Stack<SibStackData, ParticleInterface> SibStack;
-
-}
+ int GetCapacity() const { return 8000; }
+
+ void SetId(const int i, const int v) { s_plist_.llist[i] = v; }
+ void SetEnergy(const int i, const EnergyType v) { s_plist_.p[3][i] = v / 1_GeV; }
+ void SetMomentum(const int i, const super_stupid::MomentumVector& v) {
+ auto tmp = v.GetComponents();
+ for (int idx = 0; idx < 3; ++idx)
+ s_plist_.p[idx][i] = tmp[idx] / 1_GeV * constants::c;
+ }
+
+ int GetId(const int i) const { return s_plist_.llist[i]; }
+
+ EnergyType GetEnergy(const int i) const { return s_plist_.p[3][i] * 1_GeV; }
+
+ super_stupid::MomentumVector GetMomentum(const int i) const {
+ CoordinateSystem& rootCS =
+ RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+ corsika::geometry::QuantityVector<momentum_d> components{
+ s_plist_.p[0][i] * 1_GeV / constants::c,
+ s_plist_.p[1][i] * 1_GeV / constants::c,
+ s_plist_.p[2][i] * 1_GeV / constants::c};
+ super_stupid::MomentumVector v1(rootCS, components);
+ return v1;
+ }
+
+ void Copy(const int i1, const int i2) {
+ s_plist_.llist[i1] = s_plist_.llist[i2];
+ s_plist_.p[3][i1] = s_plist_.p[3][i2];
+ }
+
+ protected:
+ void IncrementSize() { s_plist_.np++; }
+ void DecrementSize() {
+ if (s_plist_.np > 0) { s_plist_.np--; }
+ }
+ };
+
+ template <typename StackIteratorInterface>
+ class ParticleInterface : public ParticleBase<StackIteratorInterface> {
+ using ParticleBase<StackIteratorInterface>::GetStackData;
+ using ParticleBase<StackIteratorInterface>::GetIndex;
+
+ public:
+ void SetEnergy(const EnergyType v) { GetStackData().SetEnergy(GetIndex(), v); }
+ EnergyType GetEnergy() const { return GetStackData().GetEnergy(GetIndex()); }
+ void SetPID(const int v) { GetStackData().SetId(GetIndex(), v); }
+ corsika::process::sibyll::SibyllCode GetPID() const {
+ return static_cast<corsika::process::sibyll::SibyllCode>(
+ GetStackData().GetId(GetIndex()));
+ }
+ super_stupid::MomentumVector GetMomentum() const {
+ return GetStackData().GetMomentum(GetIndex());
+ }
+ void SetMomentum(const super_stupid::MomentumVector& v) {
+ GetStackData().SetMomentum(GetIndex(), v);
+ }
+ };
+
+ typedef Stack<SibStackData, ParticleInterface> SibStack;
+
+} // namespace corsika::process::sibyll
#endif
diff --git a/Processes/Sibyll/sibyll2.3c.cc b/Processes/Sibyll/sibyll2.3c.cc
index 336617fbfe42e31dcc00b9d8215e8a6f5e24c485..77ad3052724c066fe474c1d41d44cb013f0786a4 100644
--- a/Processes/Sibyll/sibyll2.3c.cc
+++ b/Processes/Sibyll/sibyll2.3c.cc
@@ -2,7 +2,6 @@
#include <corsika/random/RNGManager.h>
-
double s_rndm_(int&) {
static corsika::random::RNG& rmng =
corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
diff --git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc
index 5d91a5b4b4f2ba06f93fdc25bad74c095c66fac6..a2cbb9e1b2786acab4fe950549c6a0879744aea1 100644
--- a/Processes/Sibyll/testSibyll.cc
+++ b/Processes/Sibyll/testSibyll.cc
@@ -9,8 +9,8 @@
* the license.
*/
-#include <corsika/process/sibyll/Interaction.h>
#include <corsika/process/sibyll/Decay.h>
+#include <corsika/process/sibyll/Interaction.h>
#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/particles/ParticleProperties.h>
@@ -58,12 +58,12 @@ TEST_CASE("Sibyll", "[processes]") {
REQUIRE(process::sibyll::GetSibyllXSCode(corsika::particles::Code::K0Long) == 3);
REQUIRE(process::sibyll::GetSibyllXSCode(corsika::particles::Code::SigmaPlus) == 1);
REQUIRE(process::sibyll::GetSibyllXSCode(corsika::particles::Code::PiMinus) == 2);
- }
+ }
}
#include <corsika/geometry/Point.h>
-#include <corsika/geometry/Vector.h>
#include <corsika/geometry/RootCoordinateSystem.h>
+#include <corsika/geometry/Vector.h>
#include <corsika/units/PhysicalUnits.h>
@@ -74,24 +74,26 @@ using namespace corsika::units::si;
TEST_CASE("SibyllInterface", "[processes]") {
- auto const& cs = geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
- geometry::Point const origin(cs, {0_m, 0_m, 0_m});
- geometry::Vector<corsika::units::si::SpeedType::dimension_type> v(cs, 0_m / second,
- 0_m / second, 1_m / second);
+ auto const& cs =
+ geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+ geometry::Point const origin(cs, {0_m, 0_m, 0_m});
+ geometry::Vector<corsika::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);
setup::Stack stack;
auto particle = stack.NewParticle();
-
+
SECTION("InteractionInterface") {
Interaction model;
model.Init();
- [[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(particle, stack);
- [[maybe_unused]] const double length = model.GetInteractionLength(particle, track);
-
+ [[maybe_unused]] const process::EProcessReturn ret =
+ model.DoInteraction(particle, stack);
+ [[maybe_unused]] const GrammageType length =
+ model.GetInteractionLength(particle, track);
}
SECTION("DecayInterface") {
@@ -99,10 +101,8 @@ TEST_CASE("SibyllInterface", "[processes]") {
Decay model;
model.Init();
- /*[[maybe_unused]] const process::EProcessReturn ret =*/ model.DoDecay(particle, stack);
- [[maybe_unused]] const double length = model.GetLifetime(particle);
-
+ /*[[maybe_unused]] const process::EProcessReturn ret =*/model.DoDecay(particle,
+ stack);
+ [[maybe_unused]] const TimeType time = model.GetLifetime(particle);
}
-
}
-
diff --git a/Processes/StackInspector/StackInspector.cc b/Processes/StackInspector/StackInspector.cc
index f062c6e4cf31f0a7b62853160db62b9126ace46b..52fe7e7423dc77f663a4e02c1b74c3c893ab24f8 100644
--- a/Processes/StackInspector/StackInspector.cc
+++ b/Processes/StackInspector/StackInspector.cc
@@ -44,8 +44,8 @@ process::EProcessReturn StackInspector<Stack>::DoContinuous(Particle&, setup::Tr
for (auto& iterP : s) {
EnergyType E = iterP.GetEnergy();
Etot += E;
- geometry::CoordinateSystem& rootCS =
- geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem(); // for printout
+ geometry::CoordinateSystem& rootCS = geometry::RootCoordinateSystem::GetInstance()
+ .GetRootCoordinateSystem(); // for printout
auto pos = iterP.GetPosition().GetCoordinates(rootCS);
cout << "StackInspector: i=" << setw(5) << fixed << (i++) << ", id=" << setw(30)
<< iterP.GetPID() << " E=" << setw(15) << scientific << (E / 1_GeV) << " GeV, "
@@ -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 d4b9f1faae3b2843b04acf700b6a46084b7658a6..07b987b46f8a03a2683b2d160ac86e24ac251560 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 {
@@ -32,7 +32,9 @@ namespace corsika::process {
void Init();
EProcessReturn DoContinuous(Particle&, corsika::setup::Trajectory&, Stack& s) const;
- double MaxStepLength(Particle&, corsika::setup::Trajectory&) const;
+ //~ template <typename Particle>
+ corsika::units::si::LengthType MaxStepLength(Particle&,
+ corsika::setup::Trajectory&) const;
private:
bool fReport;
diff --git a/Processes/StackInspector/testStackInspector.cc b/Processes/StackInspector/testStackInspector.cc
index 51f796a24a92439d2fc11a42c1fac07cc2f8ec30..5ff999c1d13721365f560cf5c9ba741907356651 100644
--- a/Processes/StackInspector/testStackInspector.cc
+++ b/Processes/StackInspector/testStackInspector.cc
@@ -16,8 +16,8 @@
#include <corsika/process/stack_inspector/StackInspector.h>
#include <corsika/geometry/Point.h>
-#include <corsika/geometry/Vector.h>
#include <corsika/geometry/RootCoordinateSystem.h>
+#include <corsika/geometry/Vector.h>
#include <corsika/units/PhysicalUnits.h>
@@ -28,28 +28,26 @@ using namespace corsika::units::si;
using namespace corsika::process::stack_inspector;
using namespace corsika;
-
TEST_CASE("StackInspector", "[processes]") {
- auto const& cs = geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
- geometry::Point const origin(cs, {0_m, 0_m, 0_m});
- geometry::Vector<corsika::units::si::SpeedType::dimension_type> v(cs, 0_m / second,
- 0_m / second, 1_m / second);
+ auto const& cs =
+ geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+ geometry::Point const origin(cs, {0_m, 0_m, 0_m});
+ geometry::Vector<corsika::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);
setup::Stack stack;
auto particle = stack.NewParticle();
-
+
SECTION("interface") {
StackInspector<setup::Stack> model(true);
model.Init();
- [[maybe_unused]] const process::EProcessReturn ret = model.DoContinuous(particle, track, stack);
- [[maybe_unused]] const double length = model.MaxStepLength(particle, track);
-
+ [[maybe_unused]] const process::EProcessReturn ret =
+ model.DoContinuous(particle, track, stack);
+ [[maybe_unused]] const LengthType length = model.MaxStepLength(particle, track);
}
}
-
-
diff --git a/Processes/TrackingLine/TrackingLine.h b/Processes/TrackingLine/TrackingLine.h
index 205991ecf0b2b613bc97192b15494a824636a0ea..157751ea73ab261b083857a68cd2c2444469a583 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,14 @@ 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 =