From 5500cf18a993142808e03131fba8fb48efe336c6 Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Mon, 9 Nov 2020 16:20:49 +0100
Subject: [PATCH] Major re-work of tracking.
---
CMakeLists.txt | 8 +-
Documentation/Examples/boundary_example.cc | 4 +-
Documentation/Examples/cascade_example.cc | 4 +-
.../Examples/cascade_proton_example.cc | 4 +-
Documentation/Examples/em_shower.cc | 6 +-
Documentation/Examples/hybrid_MC.cc | 78 +-
Documentation/Examples/vertical_EAS.cc | 16 +-
Environment/BaseExponential.h | 12 +-
Environment/FlatExponential.h | 4 +-
Environment/HomogeneousMedium.h | 4 +-
Environment/IMediumModel.h | 4 +-
Environment/InhomogeneousMedium.h | 4 +-
Environment/LinearApproximationIntegrator.h | 12 +-
Environment/SlidingPlanarExponential.h | 8 +-
Environment/VolumeTreeNode.h | 5 +-
Environment/testEnvironment.cc | 20 +-
Framework/Cascade/Cascade.h | 332 +--------
Framework/Cascade/testCascade.cc | 7 +-
Framework/Geometry/CMakeLists.txt | 1 +
Framework/Geometry/Helix.h | 4 +
Framework/Geometry/Intersections.hpp | 55 ++
Framework/Geometry/Line.h | 17 +-
Framework/Geometry/Sphere.h | 6 +-
Framework/Geometry/Trajectory.h | 196 ++++-
Framework/Geometry/testGeometry.cc | 11 +-
Framework/Units/PhysicalUnits.h | 4 +-
Processes/CMakeLists.txt | 3 +
.../CONEXSourceCut/testCONEXSourceCut.cc | 34 +-
.../ObservationPlane/ObservationPlane.cc | 98 +--
.../ObservationPlane/testObservationPlane.cc | 2 +-
.../StackInspector/testStackInspector.cc | 2 +-
Processes/Tracking/CMakeLists.txt | 34 +
Processes/Tracking/Intersect.hpp | 128 ++++
.../TrackingLeapFrogCurved/CMakeLists.txt | 44 ++
Processes/TrackingLeapFrogCurved/Tracking.h | 205 ++++++
.../testTrackingLeapFrogCurved.cc | 127 ++++
.../testTrackingLeapFrogStraight.cc | 127 ++++
.../TrackingLeapFrogStraight/CMakeLists.txt | 43 ++
.../Tracking.cc} | 6 +-
Processes/TrackingLeapFrogStraight/Tracking.h | 183 +++++
.../testTrackingBFieldStack.h | 58 ++
.../testTrackingLeapFrogStraight.cc | 127 ++++
Processes/TrackingLine/CMakeLists.txt | 19 +-
Processes/TrackingLine/Tracking.cc | 26 +
Processes/TrackingLine/Tracking.h | 126 ++++
Processes/TrackingLine/TrackingLine.h | 694 ------------------
.../TrackingLine/TrackingLine_interpolation.h | 193 -----
Processes/TrackingLine/dump_bh.hpp | 186 -----
Processes/TrackingLine/testTrackingLine.cc | 72 +-
.../TrackingLine/testTrackingLineStack.h | 15 +-
Setup/SetupEnvironment.h | 11 +-
Setup/SetupTrajectory.h | 36 +-
Stack/History/HistoryObservationPlane.cpp | 14 +-
do-copyright.py | 9 +-
54 files changed, 1704 insertions(+), 1744 deletions(-)
create mode 100644 Framework/Geometry/Intersections.hpp
create mode 100644 Processes/Tracking/CMakeLists.txt
create mode 100644 Processes/Tracking/Intersect.hpp
create mode 100644 Processes/TrackingLeapFrogCurved/CMakeLists.txt
create mode 100644 Processes/TrackingLeapFrogCurved/Tracking.h
create mode 100644 Processes/TrackingLeapFrogCurved/testTrackingLeapFrogCurved.cc
create mode 100644 Processes/TrackingLeapFrogCurved/testTrackingLeapFrogStraight.cc
create mode 100644 Processes/TrackingLeapFrogStraight/CMakeLists.txt
rename Processes/{TrackingLine/TrackingLine.cc => TrackingLeapFrogStraight/Tracking.cc} (92%)
create mode 100644 Processes/TrackingLeapFrogStraight/Tracking.h
create mode 100644 Processes/TrackingLeapFrogStraight/testTrackingBFieldStack.h
create mode 100644 Processes/TrackingLeapFrogStraight/testTrackingLeapFrogStraight.cc
create mode 100644 Processes/TrackingLine/Tracking.cc
create mode 100644 Processes/TrackingLine/Tracking.h
delete mode 100644 Processes/TrackingLine/TrackingLine.h
delete mode 100644 Processes/TrackingLine/TrackingLine_interpolation.h
delete mode 100644 Processes/TrackingLine/dump_bh.hpp
diff --git a/CMakeLists.txt b/CMakeLists.txt
index b126b3eed..e24a55b95 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -4,10 +4,10 @@ cmake_minimum_required (VERSION 3.9)
# prevent in-source builds and give warning message
if ("${CMAKE_BINARY_DIR}" STREQUAL "${CMAKE_SOURCE_DIR}")
message (FATAL_ERROR "In-source builds are disabled.
- Please create a subfolder and use `cmake ..` inside it.
+ Please create a build-dir and use `cmake <source-dir>` inside it.
NOTE: cmake will now create CMakeCache.txt and CMakeFiles/*.
You must delete them, or cmake will refuse to work.")
-endif ()
+endif ()
project (
corsika
@@ -47,7 +47,7 @@ set (CMAKE_CXX_STANDARD 17)
set (CMAKE_CXX_EXTENSIONS OFF)
enable_testing ()
set (CTEST_OUTPUT_ON_FAILURE 1)
-list(APPEND CMAKE_CTEST_ARGUMENTS "--output-on-failure")
+list (APPEND CMAKE_CTEST_ARGUMENTS "--output-on-failure")
# Set the possible values of build type for cmake-gui and command line check
set (ALLOWED_BUILD_TYPES Debug Release MinSizeRel RelWithDebInfo Coverage)
@@ -90,7 +90,7 @@ set (CMAKE_SHARED_LINKER_FLAGS_COVERAGE "--coverage")
add_compile_options ("$<$<OR:$<CXX_COMPILER_ID:Clang>,$<CXX_COMPILER_ID:AppleClang>>:-Wno-nonportable-include-path>")
# set a flag to inform code that we are in debug mode
-set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -DDEBUG")
+set (CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -DDEBUG")
# COAST - interface, this requires CORSIKA7 to be installed first
# COAST will allow you to program code in CORSIKA8 and execute it inside CORSIKA7
diff --git a/Documentation/Examples/boundary_example.cc b/Documentation/Examples/boundary_example.cc
index c8fa31e6e..010cb8db2 100644
--- a/Documentation/Examples/boundary_example.cc
+++ b/Documentation/Examples/boundary_example.cc
@@ -8,7 +8,7 @@
#include <corsika/cascade/Cascade.h>
#include <corsika/process/ProcessSequence.h>
-#include <corsika/process/tracking_line/TrackingLine.h>
+#include <corsika/process/tracking_line/Tracking.h>
#include <corsika/setup/SetupEnvironment.h>
#include <corsika/setup/SetupStack.h>
@@ -123,7 +123,7 @@ int main() {
universe.AddChild(std::move(world));
// setup processes, decays and interactions
- tracking_line::TrackingLine tracking;
+ tracking_line::Tracking tracking;
random::RNGManager::GetInstance().RegisterRandomStream("sibyll");
process::sibyll::Interaction sibyll;
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index f3cd9c19a..2707e3dde 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -10,7 +10,7 @@
#include <corsika/process/ProcessSequence.h>
#include <corsika/process/energy_loss/EnergyLoss.h>
#include <corsika/process/stack_inspector/StackInspector.h>
-#include <corsika/process/tracking_line/TrackingLine.h>
+#include <corsika/process/tracking_line/Tracking.h>
#include <corsika/setup/SetupEnvironment.h>
#include <corsika/setup/SetupStack.h>
@@ -135,7 +135,7 @@ int main() {
}
// setup processes, decays and interactions
- tracking_line::TrackingLine tracking;
+ tracking_line::Tracking tracking;
stack_inspector::StackInspector<setup::Stack> stackInspect(1, true, E0);
random::RNGManager::GetInstance().RegisterRandomStream("sibyll");
diff --git a/Documentation/Examples/cascade_proton_example.cc b/Documentation/Examples/cascade_proton_example.cc
index 2c5e1489e..6946f2490 100644
--- a/Documentation/Examples/cascade_proton_example.cc
+++ b/Documentation/Examples/cascade_proton_example.cc
@@ -10,7 +10,7 @@
#include <corsika/process/ProcessSequence.h>
#include <corsika/process/hadronic_elastic_model/HadronicElasticModel.h>
#include <corsika/process/stack_inspector/StackInspector.h>
-#include <corsika/process/tracking_line/TrackingLine.h>
+#include <corsika/process/tracking_line/Tracking.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
@@ -121,7 +121,7 @@ int main() {
}
// setup processes, decays and interactions
- tracking_line::TrackingLine tracking;
+ tracking_line::Tracking tracking;
stack_inspector::StackInspector<setup::Stack> stackInspect(1, true, E0);
random::RNGManager::GetInstance().RegisterRandomStream("sibyll");
diff --git a/Documentation/Examples/em_shower.cc b/Documentation/Examples/em_shower.cc
index 297ab1f67..d834afe77 100644
--- a/Documentation/Examples/em_shower.cc
+++ b/Documentation/Examples/em_shower.cc
@@ -21,7 +21,7 @@
#include <corsika/process/proposal/ContinuousProcess.h>
#include <corsika/process/proposal/Interaction.h>
#include <corsika/process/track_writer/TrackWriter.h>
-#include <corsika/process/tracking_line/TrackingLine.h>
+#include <corsika/process/tracking_line/Tracking.h>
#include <corsika/random/RNGManager.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
@@ -152,12 +152,12 @@ int main(int argc, char** argv) {
Plane const obsPlane(showerCore, Vector<dimensionless_d>(rootCS, {0., 0., 1.}));
process::observation_plane::ObservationPlane observationLevel(
- obsPlane, Vector<dimensionless_d>(rootCS, {1., 0., 0.}), "particles.dat");
+ obsPlane, Vector<dimensionless_d>(rootCS, {1., 0., 0.}), "particles.dat");
auto sequence = process::sequence(proposalCounted, em_continuous, longprof, cut,
observationLevel, trackWriter);
// define air shower object, run simulation
- tracking_line::TrackingLine tracking;
+ tracking_line::Tracking tracking;
cascade::Cascade EAS(env, tracking, sequence, stack);
// to fix the point of first interaction, uncomment the following two lines:
diff --git a/Documentation/Examples/hybrid_MC.cc b/Documentation/Examples/hybrid_MC.cc
index 33dec55b0..a9f7432ae 100644
--- a/Documentation/Examples/hybrid_MC.cc
+++ b/Documentation/Examples/hybrid_MC.cc
@@ -22,6 +22,7 @@
#include <corsika/geometry/Sphere.h>
#include <corsika/logging/Logging.h>
#include <corsika/process/ProcessSequence.h>
+#include <corsika/process/SwitchProcessSequence.h>
#include <corsika/process/StackProcess.h>
#include <corsika/process/conex_source_cut/CONEXSourceCut.h>
#include <corsika/process/energy_loss/EnergyLoss.h>
@@ -33,8 +34,7 @@
#include <corsika/process/sibyll/Decay.h>
#include <corsika/process/sibyll/Interaction.h>
#include <corsika/process/sibyll/NuclearInteraction.h>
-#include <corsika/process/switch_process/SwitchProcess.h>
-#include <corsika/process/tracking_line/TrackingLine.h>
+#include <corsika/process/tracking_line/Tracking.h>
#include <corsika/process/urqmd/UrQMD.h>
#include <corsika/random/RNGManager.h>
#include <corsika/setup/SetupStack.h>
@@ -72,14 +72,17 @@ void registerRandomStreams(const int seed) {
random::RNGManager::GetInstance().SeedAll(seed);
}
+template <typename T>
+using MyExtraEnv = environment::MediumPropertyModel<environment::UniformMagneticField<T>>;
+
int main(int argc, char** argv) {
logging::SetLevel(logging::level::info);
- C8LOG_INFO("vertical_EAS");
+ C8LOG_INFO("hybrid_MC");
if (argc < 4) {
- std::cerr << "usage: vertical_EAS <A> <Z> <energy/GeV> [seed]" << std::endl;
+ std::cerr << "usage: hybrid_MC <A> <Z> <energy/GeV> [seed]" << std::endl;
std::cerr << " if no seed is given, a random seed is chosen" << std::endl;
return 1;
}
@@ -95,37 +98,20 @@ int main(int argc, char** argv) {
EnvType env;
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
Point const center{rootCS, 0_m, 0_m, 0_m};
- environment::LayeredSphericalAtmosphereBuilder<setup::EnvironmentInterface> builder{
- center};
+ auto builder = environment::make_layered_spherical_atmosphere_builder<
+ setup::EnvironmentInterface,
+ MyExtraEnv>::create(center, units::constants::EarthRadius::Mean,
+ environment::Medium::AirDry1Atm,
+ geometry::Vector{rootCS, 0_T, 0_T, 1_T});
builder.setNuclearComposition(
{{particles::Code::Nitrogen, particles::Code::Oxygen},
{0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now
- builder.addExponentialLayer<
- environment::UniformMediumType<environment::UniformMagneticField<
- SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
- 1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km, environment::EMediumType::eAir,
- geometry::Vector(rootCS, 0_T, 0_T, 1_T));
- builder.addExponentialLayer<
- environment::UniformMediumType<environment::UniformMagneticField<
- SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
- 1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km, environment::EMediumType::eAir,
- geometry::Vector(rootCS, 0_T, 0_T, 1_T));
- builder.addExponentialLayer<
- environment::UniformMediumType<environment::UniformMagneticField<
- SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
- 1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km, environment::EMediumType::eAir,
- geometry::Vector(rootCS, 0_T, 0_T, 1_T));
- builder.addExponentialLayer<
- environment::UniformMediumType<environment::UniformMagneticField<
- SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
- 540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km, environment::EMediumType::eAir,
- geometry::Vector(rootCS, 0_T, 0_T, 1_T));
- builder.addLinearLayer<environment::UniformMediumType<
- environment::UniformMagneticField<HomogeneousMedium<setup::EnvironmentInterface>>>>(
- 1e9_cm, 112.8_km, environment::EMediumType::eAir,
- geometry::Vector(rootCS, 0_T, 0_T, 1_T));
-
+ builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km);
+ builder.addExponentialLayer(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km);
+ builder.addExponentialLayer(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km);
+ builder.addExponentialLayer(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km);
+ builder.addLinearLayer(1e9_cm, 112.8_km);
builder.assemble(env);
// setup particle stack, and add primary particle
@@ -222,7 +208,8 @@ int main(int argc, char** argv) {
process::energy_loss::EnergyLoss eLoss{showerAxis, cut.GetECut()};
corsika::process::conex_source_cut::CONEXSourceCut conex(
- center, showerAxis, t, injectionHeight, E0, particles::Code::Proton);
+ center, showerAxis, t, injectionHeight, E0,
+ particles::GetPDG(particles::Code::Proton));
process::on_shell_check::OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);
@@ -236,17 +223,26 @@ int main(int argc, char** argv) {
process::interaction_counter::InteractionCounter urqmdCounted{urqmd};
// assemble all processes into an ordered process list
-
- auto sibyllSequence = sibyllNucCounted << sibyllCounted;
- process::switch_process::SwitchProcess switchProcess(urqmdCounted, sibyllSequence,
- 55_GeV);
- auto decaySequence = decayPythia << decaySibyll;
-
- auto sequence = switchProcess << reset_particle_mass << decaySequence << eLoss << cut
- << conex << longprof << observationLevel;
+ struct EnergySwitch {
+ HEPEnergyType cutE_;
+ EnergySwitch(HEPEnergyType cutE)
+ : cutE_(cutE) {}
+ process::SwitchResult operator()(const setup::Stack::ParticleType& p) {
+ if (p.GetEnergy() < cutE_)
+ return process::SwitchResult::First;
+ else
+ return process::SwitchResult::Second;
+ }
+ };
+ auto hadronSequence =
+ process::select(urqmdCounted, process::sequence(sibyllNucCounted, sibyllCounted),
+ EnergySwitch(55_GeV));
+ auto decaySequence = process::sequence(decayPythia, decaySibyll);
+ auto sequence = process::sequence(hadronSequence, reset_particle_mass, decaySequence,
+ eLoss, cut, conex, longprof, observationLevel);
// define air shower object, run simulation
- tracking_line::TrackingLine tracking;
+ tracking_line::Tracking tracking;
cascade::Cascade EAS(env, tracking, sequence, stack);
// to fix the point of first interaction, uncomment the following two lines:
diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc
index f9b44e638..d5cb7baa0 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -33,13 +33,14 @@
#include <corsika/process/observation_plane/ObservationPlane.h>
#include <corsika/process/on_shell_check/OnShellCheck.h>
#include <corsika/process/particle_cut/ParticleCut.h>
+#include <corsika/process/track_writer/TrackWriter.h>
#include <corsika/process/proposal/ContinuousProcess.h>
#include <corsika/process/proposal/Interaction.h>
#include <corsika/process/pythia/Decay.h>
#include <corsika/process/sibyll/Decay.h>
#include <corsika/process/sibyll/Interaction.h>
#include <corsika/process/sibyll/NuclearInteraction.h>
-#include <corsika/process/tracking_line/TrackingLine.h>
+#include <corsika/process/tracking_line/Tracking.h>
#include <corsika/process/urqmd/UrQMD.h>
#include <corsika/random/RNGManager.h>
#include <corsika/setup/SetupStack.h>
@@ -109,7 +110,7 @@ int main(int argc, char** argv) {
setup::EnvironmentInterface,
MyExtraEnv>::create(center, units::constants::EarthRadius::Mean,
environment::Medium::AirDry1Atm,
- geometry::Vector{rootCS, 0_T, 0_T, 1_T});
+ geometry::Vector{rootCS, 0_T, 50_uT, 0_T});
builder.setNuclearComposition(
{{particles::Code::Nitrogen, particles::Code::Oxygen},
{0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now
@@ -147,11 +148,12 @@ int main(int argc, char** argv) {
cout << "input momentum: " << plab.GetComponents() / 1_GeV << ", norm = " << plab.norm()
<< endl;
- auto const observationHeight = 1.4_km + seaLevel;
- auto const injectionHeight = 112.75_km + seaLevel;
+ auto const observationHeight = 0_km + builder.getEarthRadius();
+ auto const injectionHeight = 112.75_km + builder.getEarthRadius();
auto const t = -observationHeight * cos(thetaRad) +
sqrt(-units::static_pow<2>(sin(thetaRad) * observationHeight) +
units::static_pow<2>(injectionHeight));
+
Point const showerCore{rootCS, 0_m, 0_m, observationHeight};
Point const injectionPos =
showerCore +
@@ -225,10 +227,6 @@ int main(int argc, char** argv) {
decaySibyll.PrintDecayConfig();
- process::particle_cut::ParticleCut cut{60_GeV, false, true};
- process::proposal::Interaction proposal(env, cut.GetECut());
- process::proposal::ContinuousProcess em_continuous(env, cut.GetECut());
- process::interaction_counter::InteractionCounter proposalCounted(proposal);
process::on_shell_check::OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);
@@ -263,7 +261,7 @@ int main(int argc, char** argv) {
proposalCounted, em_continuous, cut, observationLevel, longprof);
// define air shower object, run simulation
- tracking_line::TrackingLine tracking;
+ tracking_line::Tracking tracking;
cascade::Cascade EAS(env, tracking, sequence, stack);
// to fix the point of first interaction, uncomment the following two lines:
diff --git a/Environment/BaseExponential.h b/Environment/BaseExponential.h
index 97e2689bc..c5ac3181f 100644
--- a/Environment/BaseExponential.h
+++ b/Environment/BaseExponential.h
@@ -47,12 +47,12 @@ namespace corsika::environment {
*/
// clang-format on
units::si::GrammageType IntegratedGrammage(
- geometry::Trajectory<geometry::Line> const& vLine, units::si::LengthType vL,
+ geometry::LineTrajectory const& vLine, units::si::LengthType vL,
geometry::Vector<units::si::dimensionless_d> const& vAxis) const {
if (vL == units::si::LengthType::zero()) { return units::si::GrammageType::zero(); }
- auto const uDotA = vLine.NormalizedDirection().dot(vAxis).magnitude();
- auto const rhoStart = GetImplementation().GetMassDensity(vLine.GetR0());
+ auto const uDotA = vLine.GetDirection(0).dot(vAxis).magnitude();
+ auto const rhoStart = GetImplementation().GetMassDensity(vLine.GetLine().GetR0());
if (uDotA == 0) {
return vL * rhoStart;
@@ -80,11 +80,11 @@ namespace corsika::environment {
*/
// clang-format on
units::si::LengthType ArclengthFromGrammage(
- geometry::Trajectory<geometry::Line> const& vLine,
+ geometry::LineTrajectory const& vLine,
units::si::GrammageType vGrammage,
geometry::Vector<units::si::dimensionless_d> const& vAxis) const {
- auto const uDotA = vLine.NormalizedDirection().dot(vAxis).magnitude();
- auto const rhoStart = GetImplementation().GetMassDensity(vLine.GetR0());
+ auto const uDotA = vLine.GetDirection(0).dot(vAxis).magnitude();
+ auto const rhoStart = GetImplementation().GetMassDensity(vLine.GetLine().GetR0());
if (uDotA == 0) {
return vGrammage / rhoStart;
diff --git a/Environment/FlatExponential.h b/Environment/FlatExponential.h
index e1b0cd5b6..c29dd91c7 100644
--- a/Environment/FlatExponential.h
+++ b/Environment/FlatExponential.h
@@ -52,13 +52,13 @@ namespace corsika::environment {
NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
units::si::GrammageType IntegratedGrammage(
- geometry::Trajectory<geometry::Line> const& vLine,
+ geometry::LineTrajectory const& vLine,
units::si::LengthType vTo) const override {
return Base::IntegratedGrammage(vLine, vTo, fAxis);
}
units::si::LengthType ArclengthFromGrammage(
- geometry::Trajectory<geometry::Line> const& vLine,
+ geometry::LineTrajectory const& vLine,
units::si::GrammageType vGrammage) const override {
return Base::ArclengthFromGrammage(vLine, vGrammage, fAxis);
}
diff --git a/Environment/HomogeneousMedium.h b/Environment/HomogeneousMedium.h
index 49f4cb80d..d630d8a6b 100644
--- a/Environment/HomogeneousMedium.h
+++ b/Environment/HomogeneousMedium.h
@@ -42,14 +42,14 @@ namespace corsika::environment {
NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
corsika::units::si::GrammageType IntegratedGrammage(
- corsika::geometry::Trajectory<corsika::geometry::Line> const&,
+ corsika::geometry::LineTrajectory const&,
corsika::units::si::LengthType pTo) const override {
using namespace corsika::units::si;
return pTo * fDensity;
}
corsika::units::si::LengthType ArclengthFromGrammage(
- corsika::geometry::Trajectory<corsika::geometry::Line> const&,
+ corsika::geometry::LineTrajectory const&,
corsika::units::si::GrammageType pGrammage) const override {
return pGrammage / fDensity;
}
diff --git a/Environment/IMediumModel.h b/Environment/IMediumModel.h
index da6c9ca47..388f79c43 100644
--- a/Environment/IMediumModel.h
+++ b/Environment/IMediumModel.h
@@ -26,11 +26,11 @@ namespace corsika::environment {
// todo: think about the mixin inheritance of the trajectory vs the BaseTrajectory
// approach; for now, only lines are supported
virtual corsika::units::si::GrammageType IntegratedGrammage(
- corsika::geometry::Trajectory<corsika::geometry::Line> const&,
+ corsika::geometry::LineTrajectory const&,
corsika::units::si::LengthType) const = 0;
virtual corsika::units::si::LengthType ArclengthFromGrammage(
- corsika::geometry::Trajectory<corsika::geometry::Line> const&,
+ corsika::geometry::LineTrajectory const&,
corsika::units::si::GrammageType) const = 0;
virtual NuclearComposition const& GetNuclearComposition() const = 0;
diff --git a/Environment/InhomogeneousMedium.h b/Environment/InhomogeneousMedium.h
index 60f591923..87002d6b1 100644
--- a/Environment/InhomogeneousMedium.h
+++ b/Environment/InhomogeneousMedium.h
@@ -41,13 +41,13 @@ namespace corsika::environment {
NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
corsika::units::si::GrammageType IntegratedGrammage(
- corsika::geometry::Trajectory<corsika::geometry::Line> const& pLine,
+ corsika::geometry::LineTrajectory const& pLine,
corsika::units::si::LengthType pTo) const override {
return fDensityFunction.IntegrateGrammage(pLine, pTo);
}
corsika::units::si::LengthType ArclengthFromGrammage(
- corsika::geometry::Trajectory<corsika::geometry::Line> const& pLine,
+ corsika::geometry::LineTrajectory const& pLine,
corsika::units::si::GrammageType pGrammage) const override {
return fDensityFunction.ArclengthFromGrammage(pLine, pGrammage);
}
diff --git a/Environment/LinearApproximationIntegrator.h b/Environment/LinearApproximationIntegrator.h
index 1224181f3..66bef41e2 100644
--- a/Environment/LinearApproximationIntegrator.h
+++ b/Environment/LinearApproximationIntegrator.h
@@ -21,29 +21,29 @@ namespace corsika::environment {
public:
auto IntegrateGrammage(
- corsika::geometry::Trajectory<corsika::geometry::Line> const& line,
+ corsika::geometry::LineTrajectory const& line,
corsika::units::si::LengthType length) const {
auto const c0 = GetImplementation().EvaluateAt(line.GetPosition(0));
auto const c1 = GetImplementation().fRho.FirstDerivative(
- line.GetPosition(0), line.NormalizedDirection());
+ line.GetPosition(0), line.GetDirection(0));
return (c0 + 0.5 * c1 * length) * length;
}
auto ArclengthFromGrammage(
- corsika::geometry::Trajectory<corsika::geometry::Line> const& line,
+ corsika::geometry::LineTrajectory const& line,
corsika::units::si::GrammageType grammage) const {
auto const c0 = GetImplementation().fRho(line.GetPosition(0));
auto const c1 = GetImplementation().fRho.FirstDerivative(
- line.GetPosition(0), line.NormalizedDirection());
+ line.GetPosition(0), line.GetDirection(0));
return (1 - 0.5 * grammage * c1 / (c0 * c0)) * grammage / c0;
}
- auto MaximumLength(corsika::geometry::Trajectory<corsika::geometry::Line> const& line,
+ auto MaximumLength(corsika::geometry::LineTrajectory const& line,
[[maybe_unused]] double relError) const {
using namespace corsika::units::si;
[[maybe_unused]] auto const c1 = GetImplementation().fRho.SecondDerivative(
- line.GetPosition(0), line.NormalizedDirection());
+ line.GetPosition(0), line.GetDirection(0));
// todo: provide a real, working implementation
return 1_m * std::numeric_limits<double>::infinity();
diff --git a/Environment/SlidingPlanarExponential.h b/Environment/SlidingPlanarExponential.h
index 41058d864..b37fff2d1 100644
--- a/Environment/SlidingPlanarExponential.h
+++ b/Environment/SlidingPlanarExponential.h
@@ -55,16 +55,16 @@ namespace corsika::environment {
NuclearComposition const& GetNuclearComposition() const override { return nuclComp_; }
units::si::GrammageType IntegratedGrammage(
- geometry::Trajectory<geometry::Line> const& line,
+ geometry::LineTrajectory const& line,
units::si::LengthType l) const override {
- auto const axis = (line.GetR0() - Base::fP0).normalized();
+ auto const axis = (line.GetLine().GetR0() - Base::fP0).normalized();
return Base::IntegratedGrammage(line, l, axis);
}
units::si::LengthType ArclengthFromGrammage(
- geometry::Trajectory<geometry::Line> const& line,
+ geometry::LineTrajectory const& line,
units::si::GrammageType grammage) const override {
- auto const axis = (line.GetR0() - Base::fP0).normalized();
+ auto const axis = (line.GetLine().GetR0() - Base::fP0).normalized();
return Base::ArclengthFromGrammage(line, grammage, axis);
}
};
diff --git a/Environment/VolumeTreeNode.h b/Environment/VolumeTreeNode.h
index 3acb988e8..1cd5cf222 100644
--- a/Environment/VolumeTreeNode.h
+++ b/Environment/VolumeTreeNode.h
@@ -21,6 +21,7 @@ namespace corsika::environment {
class VolumeTreeNode {
public:
using IModelProperties = TModelProperties;
+ using VTN_type = VolumeTreeNode<IModelProperties>;
using VTNUPtr = std::unique_ptr<VolumeTreeNode<IModelProperties>>;
using IMPSharedPtr = std::shared_ptr<IModelProperties>;
using VolUPtr = std::unique_ptr<corsika::geometry::Volume>;
@@ -91,8 +92,8 @@ namespace corsika::environment {
void ExcludeOverlapWith(VTNUPtr const& pNode) {
fExcludedNodes.push_back(pNode.get());
}
-
- auto* GetParent() const { return fParentNode; };
+
+ const VTN_type* GetParent() const { return fParentNode; };
auto const& GetChildNodes() const { return fChildNodes; }
diff --git a/Environment/testEnvironment.cc b/Environment/testEnvironment.cc
index e93eb64f6..37ba3049a 100644
--- a/Environment/testEnvironment.cc
+++ b/Environment/testEnvironment.cc
@@ -61,7 +61,7 @@ TEST_CASE("FlatExponential") {
SECTION("horizontal") {
Line const line(gOrigin, Vector<SpeedType::dimension_type>(
gCS, {20_cm / second, 0_m / second, 0_m / second}));
- Trajectory<Line> const trajectory(line, tEnd);
+ LineTrajectory const trajectory(line, tEnd);
CHECK((medium.IntegratedGrammage(trajectory, 2_m) / (rho0 * 2_m)) == Approx(1));
CHECK((medium.ArclengthFromGrammage(trajectory, rho0 * 5_m) / 5_m) == Approx(1));
@@ -70,7 +70,7 @@ TEST_CASE("FlatExponential") {
SECTION("vertical") {
Line const line(gOrigin, Vector<SpeedType::dimension_type>(
gCS, {0_m / second, 0_m / second, 5_m / second}));
- Trajectory<Line> const trajectory(line, tEnd);
+ LineTrajectory const trajectory(line, tEnd);
LengthType const length = 2 * lambda;
GrammageType const exact = rho0 * lambda * (exp(length / lambda) - 1);
@@ -81,11 +81,11 @@ TEST_CASE("FlatExponential") {
SECTION("escape grammage") {
Line const line(gOrigin, Vector<SpeedType::dimension_type>(
gCS, {0_m / second, 0_m / second, -5_m / second}));
- Trajectory<Line> const trajectory(line, tEnd);
+ LineTrajectory const trajectory(line, tEnd);
GrammageType const escapeGrammage = rho0 * lambda;
- CHECK(trajectory.NormalizedDirection().dot(axis).magnitude() < 0);
+ CHECK(trajectory.GetDirection(0).dot(axis).magnitude() < 0);
CHECK(medium.ArclengthFromGrammage(trajectory, 1.2 * escapeGrammage) ==
std::numeric_limits<typename GrammageType::value_type>::infinity() * 1_m);
}
@@ -93,7 +93,7 @@ TEST_CASE("FlatExponential") {
SECTION("inclined") {
Line const line(gOrigin, Vector<SpeedType::dimension_type>(
gCS, {0_m / second, 5_m / second, 5_m / second}));
- Trajectory<Line> const trajectory(line, tEnd);
+ LineTrajectory const trajectory(line, tEnd);
double const cosTheta = M_SQRT1_2;
LengthType const length = 2 * lambda;
GrammageType const exact =
@@ -127,7 +127,7 @@ TEST_CASE("SlidingPlanarExponential") {
Line const line({gCS, {0_m, 0_m, 1_m}},
Vector<SpeedType::dimension_type>(
gCS, {0_m / second, 0_m / second, 5_m / second}));
- Trajectory<Line> const trajectory(line, tEnd);
+ LineTrajectory const trajectory(line, tEnd);
CHECK(medium.GetMassDensity({gCS, {0_mm, 0_m, 3_m}}).magnitude() ==
flat.GetMassDensity({gCS, {0_mm, 0_m, 3_m}}).magnitude());
@@ -166,7 +166,7 @@ TEST_CASE("InhomogeneousMedium") {
gCS, {20_m / second, 0_m / second, 0_m / second}));
auto const tEnd = 5_s;
- Trajectory<Line> const trajectory(line, tEnd);
+ LineTrajectory const trajectory(line, tEnd);
Exponential const e;
DensityFunction<decltype(e), LinearApproximationIntegrator> const rho(e);
@@ -292,7 +292,7 @@ TEST_CASE("UniformMagneticField w/ Homogeneous Medium") {
auto const tEnd = 1_s;
// and the associated trajectory
- Trajectory<Line> const trajectory(line, tEnd);
+ LineTrajectory const trajectory(line, tEnd);
// and check the integrated grammage
CHECK((medium.IntegratedGrammage(trajectory, 3_m) / (density * 3_m)) == Approx(1));
@@ -395,7 +395,7 @@ TEST_CASE("UniformRefractiveIndex w/ Homogeneous") {
auto const tEnd = 1_s;
// and the associated trajectory
- Trajectory<Line> const trajectory(line, tEnd);
+ LineTrajectory const trajectory(line, tEnd);
// and check the integrated grammage
CHECK((medium.IntegratedGrammage(trajectory, 3_m) / (density * 3_m)) == Approx(1));
@@ -469,7 +469,7 @@ TEST_CASE("MediumPropertyModel w/ Homogeneous") {
auto const tEnd = 1_s;
// and the associated trajectory
- Trajectory<Line> const trajectory(line, tEnd);
+ LineTrajectory const trajectory(line, tEnd);
// and check the integrated grammage
CHECK((medium.IntegratedGrammage(trajectory, 3_m) / (density * 3_m)) == Approx(1));
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index a64298060..9ae6724df 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -37,53 +37,6 @@
using boost::typeindex::type_id_with_cvr;
#include <fstream>
-#include <boost/histogram.hpp>
-#include <boost/histogram/ostream.hpp>
-#include <corsika/process/tracking_line/dump_bh.hpp>
-using namespace boost::histogram;
-/*static auto histL2 = make_histogram(axis::regular<>(100, 0, 60000, "L'"));
-static auto histS2 = make_histogram(axis::regular<>(100, 0, 60000, "S"));
-static auto histB2 = make_histogram(axis::regular<>(100, 0, 60000, "Bogenlänge"));*/
-static auto histLlog2 = make_histogram(axis::regular<double, axis::transform::log>(100, 1e-3, 1e7, "Leap-Frog-length L'"));
-static auto histLlog2int = make_histogram(axis::regular<double, axis::transform::log>(100, 1e-3, 1e7, "Leap-Frog-length L'"));
-static auto histLlog2dec = make_histogram(axis::regular<double, axis::transform::log>(100, 1e-3, 1e7, "Leap-Frog-length L'"));
-static auto histLlog2max = make_histogram(axis::regular<double, axis::transform::log>(100, 1e-3, 1e7, "Leap-Frog-length L'"));
-static auto histLlog2geo = make_histogram(axis::regular<double, axis::transform::log>(100, 1e-3, 1e7, "Leap-Frog-length L'"));
-static auto histLlog2mag = make_histogram(axis::regular<double, axis::transform::log>(100, 1e-3, 1e7, "Leap-Frog-length L'"));
-
-/*static auto histSlog2 = make_histogram(axis::regular<double, axis::transform::log>(100, 1e-3, 1e7, "Direct Length S"));
-static auto histBlog2 = make_histogram(axis::regular<double, axis::transform::log>(100, 1e-3, 1e7, "Arc Length B"));
-static auto histLB2 = make_histogram(axis::regular<>(100, 0, 0.01, "L - B"));
-static auto histLS2 = make_histogram(axis::regular<>(100, 0, 0.01, "L - S"));
-static auto histLBrel2 = make_histogram(axis::regular<double, axis::transform::log> (20,1e-11,1e-6,"L/B -1"));
-static auto histLSrel2 = make_histogram(axis::regular<double, axis::transform::log>(20,1e-11,1e-6, "L/S -1"));
-static auto histELSrel2 = make_histogram(axis::regular<double, axis::transform::log>(20,1e-11,1e-6, "L/S -1"),axis::regular<double, axis::transform::log>(20, 0.1, 1e4, "E / GeV"));
-static auto histBS2 = make_histogram(axis::regular<>(100, 0, 0.01, "B - S")); */
-
-static auto histLp2 = make_histogram(axis::regular<double, axis::transform::log>(100, 1e-3, 1e7, "L' für Protonen"));
-static auto histLpi2 = make_histogram(axis::regular<double, axis::transform::log>(100, 1e-3, 1e7, "L' für Pionen"));
-static auto histLpi2int = make_histogram(axis::regular<double, axis::transform::log>(100, 1e-3, 1e7, "Leap-Frog-length L'"));
-static auto histLpi2dec = make_histogram(axis::regular<double, axis::transform::log>(100, 1e-3, 1e7, "Leap-Frog-length L'"));
-static auto histLpi2max = make_histogram(axis::regular<double, axis::transform::log>(100, 1e-3, 1e7, "Leap-Frog-length L'"));
-static auto histLpi2geo = make_histogram(axis::regular<double, axis::transform::log>(100, 1e-3, 1e7, "Leap-Frog-length L'"));
-static auto histLpi2mag = make_histogram(axis::regular<double, axis::transform::log>(100, 1e-3, 1e7, "Leap-Frog-length L'"));
-static auto histLmu2 = make_histogram(axis::regular<double, axis::transform::log>(100, 1e-3, 1e7, "L' für Myonen"));
-static auto histLmu2int = make_histogram(axis::regular<double, axis::transform::log>(100, 1e-3, 1e7, "Leap-Frog-length L'"));
-static auto histLmu2dec = make_histogram(axis::regular<double, axis::transform::log>(100, 1e-3, 1e7, "Leap-Frog-length L'"));
-static auto histLmu2max = make_histogram(axis::regular<double, axis::transform::log>(100, 1e-3, 1e7, "Leap-Frog-length L'"));
-static auto histLmu2geo = make_histogram(axis::regular<double, axis::transform::log>(100, 1e-3, 1e7, "Leap-Frog-length L'"));
-static auto histLmu2mag = make_histogram(axis::regular<double, axis::transform::log>(100, 1e-3, 1e7, "Leap-Frog-length L'"));
-static auto histLe2 = make_histogram(axis::regular<double, axis::transform::log>(100, 1e-3, 1e7, "L' für Elektronen"));
-static auto histLy2 = make_histogram(axis::regular<double, axis::transform::log>(100, 1e-3, 1e7, "L' für Photonen"));
-
-static double stepradius = 0;
-static int N = 0;
-static double stepradiusp = 0;
-static int Np = 0;
-static double stepradiuspi = 0;
-static int Npi = 0;
-static double stepradiusmu = 0;
-static int Nmu = 0;
/**
@@ -125,17 +78,6 @@ namespace corsika::cascade {
std::remove_pointer_t<decltype(((Particle*)nullptr)->GetNode())>;
using MediumInterface = typename VolumeTreeNode::IModelProperties;
- private:
- // Data members
- corsika::environment::Environment<MediumInterface> const& environment_;
- TTracking& tracking_;
- TProcessList& process_sequence_;
- TStack& stack_;
- corsika::random::RNG& rng_ =
- corsika::random::RNGManager::GetInstance().GetRandomStream("cascade");
- unsigned int count_ = 0;
-
- private:
// we only want fully configured objects
Cascade() = delete;
@@ -158,132 +100,8 @@ namespace corsika::cascade {
}
~Cascade(){
- /*std::ofstream myfile;
- myfile.open ("stepradius.txt");
- myfile << "All charged particles " << stepradius/N << std::endl;
- myfile << "Protons " << stepradiusp/Np << std::endl;
- myfile << "Pions " << stepradiuspi/Npi << std::endl;
- myfile << "Muons " << stepradiusmu/Nmu << std::endl;
- myfile.close();
-
- /*std::cout << histLBrel << std::endl;
- std::cout << histLSrel << std::endl;*/
-
-
- /*std::ofstream file1("histL2.json");
- dump_bh(file1, histL2);
- file1.close();
- std::ofstream file2("histS2.json");
- dump_bh(file2, histS2);
- file2.close();
- std::ofstream file3("histB2.json");
- dump_bh(file3, histB2);
- file3.close();
- std::ofstream file4("histLB.json");
- dump_bh(file4, histLB);
- file4.close();
- std::ofstream file5("histLS.json");
- dump_bh(file5, histLS);
- file5.close();
- std::ofstream file6("histBS.json");
- dump_bh(file6, histBS);
- file6.close();
- std::ofstream file7("histLBrel.json");
- dump_bh(file7, histLBrel);
- file7.close();
- std::ofstream file8("histLSrel.json");
- dump_bh(file8, histLSrel);
- file8.close();
-
- std::ofstream file10("histELSrel.json");
- dump_bh(file10, histELSrel);
- file10.close(); //hier ende
- std::ofstream file19("histLy2.json");
- dump_bh(file19, histLy2);
- file19.close();
- std::ofstream file10("histLe2.json");
- dump_bh(file10, histLe2);
- file10.close();
- std::ofstream file11("histLmu2.json");
- dump_bh(file11, histLmu2);
- file11.close();
- std::ofstream file12("histLpi2.json");
- dump_bh(file12, histLpi2);
- file12.close();
- std::ofstream file13("histLp2.json");
- dump_bh(file13, histLp2);
- file13.close();
- std::ofstream file14("histLlog2.json");
- dump_bh(file14, histLlog2);
- file14.close();
- /*std::ofstream file15("histBlog2.json");
- dump_bh(file15, histBlog2);
- file15.close();
- std::ofstream file16("histSlog2.json");
- dump_bh(file16, histSlog2);
- file16.close(); //hier ende
- std::ofstream file17("histLlog2int.json");
- dump_bh(file17, histLlog2int);
- file17.close();
- std::ofstream file18("histLlog2dec.json");
- dump_bh(file18, histLlog2dec);
- file18.close();
-
- std::ofstream file20("histLlog2mag.json");
- dump_bh(file20, histLlog2mag);
- file20.close();
- std::ofstream file21("histLlog2geo.json");
- dump_bh(file21, histLlog2geo);
- file21.close();
- std::ofstream file22("histLlog2max.json");
- dump_bh(file22, histLlog2max);
- file22.close();
-
- std::ofstream filepi1("histLpi2int.json");
- dump_bh(filepi1, histLpi2int);
- filepi1.close();
- std::ofstream filepi2("histLpi2dec.json");
- dump_bh(filepi2, histLpi2dec);
- filepi2.close();
- std::ofstream filepi3("histLpi2mag.json");
- dump_bh(filepi3, histLpi2mag);
- filepi3.close();
- std::ofstream filepi4("histLpi2geo.json");
- dump_bh(filepi4, histLpi2geo);
- filepi4.close();
- std::ofstream filepi5("histLpi2max.json");
- dump_bh(filepi5, histLpi2max);
- filepi5.close();
-
- std::ofstream filemu1("histLmu2int.json");
- dump_bh(filemu1, histLmu2int);
- filemu1.close();
- std::ofstream filemu2("histLmu2dec.json");
- dump_bh(filemu2, histLmu2dec);
- filemu2.close();
- std::ofstream filemu3("histLmu2mag.json");
- dump_bh(filemu3, histLmu2mag);
- filemu3.close();
- std::ofstream filemu4("histLmu2geo.json");
- dump_bh(filemu4, histLmu2geo);
- filemu4.close();
- std::ofstream filemu5("histLmu2max.json");
- dump_bh(filemu5, histLmu2max);
- filemu5.close();*/
-
};
- /**
- * set the nodes for all particles on the stack according to their numerical
- * position
- */
- void SetNodes() {
- std::for_each(fStack.begin(), fStack.end(), [&](auto& p) {
- auto const* numericalNode =
- fEnvironment.GetUniverse()->GetContainingNode(p.GetPosition());
- p.SetNode(numericalNode);
- });
- }
/**
* The Run function is the main simulation loop, which processes
@@ -341,10 +159,6 @@ namespace corsika::cascade {
using namespace corsika;
using namespace corsika::units::si;
- // determine geometric tracking
- auto [step, geomMaxLength, nextVol] = tracking_.GetTrack(vParticle);
- [[maybe_unused]] auto const& dummy_nextVol = nextVol;
-
// determine combined total interaction length (inverse)
InverseGrammageType const total_inv_lambda =
process_sequence_.GetInverseInteractionLength(vParticle);
@@ -383,8 +197,8 @@ namespace corsika::cascade {
vParticle.GetEnergy() * units::constants::c;
// determine geometric tracking
- auto [step, geomMaxLength, magMaxLength, nextVol] = fTracking.GetTrack(vParticle);
- [[maybe_unused]] auto const& dummy_nextVol = nextVol;
+ auto [step, nextVol] = tracking_.GetTrack(vParticle);
+ auto geomMaxLength = step.GetLength(1);
// convert next_step from grammage to length
LengthType const distance_interact =
@@ -392,138 +206,28 @@ namespace corsika::cascade {
next_interact);
// determine the maximum geometric step length
- LengthType const distance_max = fProcessSequence.MaxStepLength(vParticle, stepWithoutB);
+ LengthType const distance_max = process_sequence_.MaxStepLength(vParticle, step);
C8LOG_DEBUG("distance_max={} m", distance_max / 1_m);
// take minimum of geometry, interaction, decay for next step
auto min_distance = std::min(
- {distance_interact, distance_decay, distance_max, geomMaxLength, magMaxLength});
+ {distance_interact, distance_decay, distance_max, geomMaxLength});
C8LOG_DEBUG("transport particle by : {} m "
- "Max Displacement after: {} m "
"Medium transition after: {} m "
"Decay after: {} m "
"Interaction after: {} m",
- min_distance/1_m, magMaxLength/1_m, geomMaxLength/1_m, distance_decay/1_m, distance_interact/1_m);
+ min_distance/1_m, geomMaxLength/1_m, distance_decay/1_m, distance_interact/1_m);
- // determine steplength for the magnetic field
- // because Steplength should not be min_distance
-
- auto [position, direction, L2] = fTracking.MagneticStep(vParticle, min_distance);
-
-
- int chargeNumber;
- if (corsika::particles::IsNucleus(vParticle.GetPID())) {
- chargeNumber = vParticle.GetNuclearZ();
- } else {
- chargeNumber = corsika::particles::GetChargeNumber(vParticle.GetPID());
- }
- //histL2(L2);
- histLlog2(L2);
- int pdg = static_cast<int>(particles::GetPDG(vParticle.GetPID()));
- if (min_distance == distance_interact){
- histLlog2int(L2);
- if (abs(pdg) == 13)
- histLmu2int(L2);
-
- if (abs(pdg) == 211 || abs(pdg) == 111)
- histLpi2int(L2);
- }
- if (min_distance == distance_decay) {
- histLlog2dec(L2);
- if (abs(pdg) == 13)
- histLmu2dec(L2);
- if (abs(pdg) == 211 || abs(pdg) == 111)
- histLpi2dec(L2);
- }
- if (min_distance == distance_max) {
- histLlog2max(L2);
- if (abs(pdg) == 13)
- histLmu2max(L2);
- if (abs(pdg) == 211 || abs(pdg) == 111)
- histLpi2max(L2);
- }
- if (min_distance == geomMaxLength) {
- histLlog2geo(L2);
- if (abs(pdg) == 13)
- histLmu2geo(L2);
- if (abs(pdg) == 211 || abs(pdg) == 111)
- histLpi2geo(L2);
- }
- if (min_distance == magMaxLength) {
- histLlog2mag(L2);
- if (abs(pdg) == 13)
- histLmu2mag(L2);
- if (abs(pdg) == 211 || abs(pdg) == 111)
- histLpi2mag(L2);
- }
- if (abs(pdg) == 13)
- histLmu2(L2);
- if (abs(pdg) == 11)
- histLe2(L2);
- if (abs(pdg) == 22)
- histLy2(L2);
- if (abs(pdg) == 211 || abs(pdg) == 111)
- histLpi2(L2);
- if (abs(pdg) == 2212 || abs(pdg) == 2112)
- histLp2(L2);
-
-
-
- //int chargeNumber = 0;
- if (corsika::particles::IsNucleus(vParticle.GetPID())) {
- chargeNumber = vParticle.GetNuclearZ();
- } else {
- chargeNumber = corsika::particles::GetChargeNumber(vParticle.GetPID());
- }
- if(chargeNumber != 0) {
- auto const* currentLogicalVolumeNode = vParticle.GetNode();
- auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(vParticle.GetPosition());
- geometry::Vector<SpeedType::dimension_type> velocity =
- vParticle.GetMomentum() / vParticle.GetEnergy() * corsika::units::constants::c;
- geometry::Vector<SpeedType::dimension_type> const velocityVerticalMag = velocity -
- velocity.parallelProjectionOnto(magneticfield);
- LengthType const gyroradius = vParticle.GetEnergy() * velocityVerticalMag.GetNorm() * 1_V /
- (corsika::units::constants::cSquared * abs(chargeNumber) *
- magneticfield.GetNorm() * 1_eV);
- stepradius = stepradius + min_distance/gyroradius;
- N ++;
- if (abs(pdg) == 13) {
- stepradiusmu += min_distance/gyroradius;
- Nmu ++;
- }
- if (abs(pdg) == 211 || abs(pdg) == 111) {
- stepradiuspi += min_distance/gyroradius;
- Npi ++;
- }
- if (abs(pdg) == 2212 || abs(pdg) == 2112) {
- stepradiusp += min_distance/gyroradius;
- Np ++;
- }
- }
-
-
- auto distance = position - vParticle.GetPosition();
-
- //Building Trajectory for Continuous processes
- //could also be done in MagneticStep
- geometry::Vector<SpeedType::dimension_type> velocity =
- vParticle.GetMomentum() / vParticle.GetEnergy() * corsika::units::constants::c;
- if (distance.norm() != 0_m) {
- velocity = distance.normalized() * velocity.norm();
- }
- geometry::Line line(vParticle.GetPosition(), velocity);
- geometry::Trajectory<geometry::Line> stepNew(line, distance.norm() / line.GetV0().norm());
-
// here the particle is actually moved along the trajectory to new position:
- // std::visit(setup::ParticleUpdate<Particle>{vParticle}, step);
- vParticle.SetMomentum(direction * vParticle.GetMomentum().norm());
- vParticle.SetPosition(position);
- vParticle.SetTime(vParticle.GetTime() + distance.norm() / velocity.norm());
+ step.SetLength(min_distance);
+ vParticle.SetPosition(step.GetPosition(1));
+ vParticle.SetMomentum(step.GetDirection(1)*vParticle.GetMomentum().norm());
+ vParticle.SetTime(vParticle.GetTime() + step.GetDuration());
std::cout << "New Position: " << vParticle.GetPosition().GetCoordinates() << std::endl;
// apply all continuous processes on particle + track
- process::EProcessReturn status = fProcessSequence.DoContinuous(vParticle, stepNew);
+ process::EProcessReturn status = process_sequence_.DoContinuous(vParticle, step);
if (status == process::EProcessReturn::eParticleAbsorbed) {
C8LOG_DEBUG("Cascade: delete absorbed particle PID={} E={} GeV",
@@ -544,7 +248,7 @@ namespace corsika::cascade {
TStackView secondaries(vParticle);
- if (min_distance != distance_max && min_distance != magMaxLength) {
+ if (min_distance != distance_max) {
/*
Create SecondaryView object on Stack. The data container
remains untouched and identical, and 'projectil' is identical
@@ -665,6 +369,18 @@ Y8, Y8, ,8P 88 `8b `8b 88 88P Y8b d8"
Y8a. .a8P Y8a. .a8P 88 `8b Y8a a8P 88 88 "88, d8' `8b Y8a a8P
`"Y8888Y"' `"Y8888Y"' 88 `8b "Y88888P" 88 88 Y8b d8' `8b "Y88888P"
)V0G0N";
- };
+
+ private:
+ // Data members
+ corsika::environment::Environment<MediumInterface> const& environment_;
+ TTracking& tracking_;
+ TProcessList& process_sequence_;
+ TStack& stack_;
+ corsika::random::RNG& rng_ =
+ corsika::random::RNGManager::GetInstance().GetRandomStream("cascade");
+ unsigned int count_ = 0;
+
+
+ }; // end class Cascade
} // namespace corsika::cascade
diff --git a/Framework/Cascade/testCascade.cc b/Framework/Cascade/testCascade.cc
index f0514598e..1018bb002 100644
--- a/Framework/Cascade/testCascade.cc
+++ b/Framework/Cascade/testCascade.cc
@@ -13,7 +13,8 @@
#include <corsika/process/ProcessSequence.h>
#include <corsika/process/NullModel.h>
#include <corsika/process/stack_inspector/StackInspector.h>
-#include <corsika/process/tracking_line/TrackingLine.h>
+#include <corsika/process/tracking_line/Tracking.h>
+//#include <corsika/process/tracking_bfield/Tracking.h>
#include <corsika/particles/ParticleProperties.h>
@@ -132,7 +133,7 @@ TEST_CASE("Cascade", "[Cascade]") {
auto env = MakeDummyEnv();
auto const& rootCS = env.GetCoordinateSystem();
- tracking_line::TrackingLine tracking;
+ tracking_line::Tracking tracking;
stack_inspector::StackInspector<TestCascadeStack> stackInspect(1, true, E0);
process::NullModel nullModel;
@@ -152,7 +153,7 @@ TEST_CASE("Cascade", "[Cascade]") {
particles::GetMass(particles::Code::Electron)))}),
Point(rootCS, {0_m, 0_m, 10_km}), 0_ns));
- cascade::Cascade<tracking_line::TrackingLine, decltype(sequence), TestCascadeStack,
+ cascade::Cascade<tracking_line::Tracking, decltype(sequence), TestCascadeStack,
TestCascadeStackView>
EAS(env, tracking, sequence, stack);
diff --git a/Framework/Geometry/CMakeLists.txt b/Framework/Geometry/CMakeLists.txt
index de2f6c5bb..d302289b0 100644
--- a/Framework/Geometry/CMakeLists.txt
+++ b/Framework/Geometry/CMakeLists.txt
@@ -18,6 +18,7 @@ set (
QuantityVector.h
Trajectory.h
FourVector.h
+ Intersections.hpp
)
set (
diff --git a/Framework/Geometry/Helix.h b/Framework/Geometry/Helix.h
index 01fd0e0f2..5deaa0822 100644
--- a/Framework/Geometry/Helix.h
+++ b/Framework/Geometry/Helix.h
@@ -52,6 +52,10 @@ namespace corsika::geometry {
(vPerp * (cos(omegaC * t) - 1) + uPerp * sin(omegaC * t)) / omegaC;
}
+ VelocityVec GetVelocity(corsika::units::si::TimeType t) const {
+ return vPar + (vPerp * (cos(omegaC * t) - 1) + uPerp * sin(omegaC * t));
+ }
+
Point PositionFromArclength(corsika::units::si::LengthType l) const {
return GetPosition(TimeFromArclength(l));
}
diff --git a/Framework/Geometry/Intersections.hpp b/Framework/Geometry/Intersections.hpp
new file mode 100644
index 000000000..1d0ca566d
--- /dev/null
+++ b/Framework/Geometry/Intersections.hpp
@@ -0,0 +1,55 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#pragma once
+
+#include <corsika/units/PhysicalUnits.h>
+
+#include <map> // for pair
+
+namespace corsika::geometry {
+
+ /**
+ * \class Intersection
+ *
+ * Container to store and return a list of intersections of a
+ * trajectory with a geometric volume objects in space.
+ *
+ **/
+
+ class Intersections {
+
+ Intersections(const Intersections&) = delete;
+ Intersections(Intersections&&) = delete;
+ Intersections& operator=(const Intersections&) = delete;
+
+ public:
+ Intersections()
+ : has_intersections_(false) {}
+ Intersections(corsika::units::si::TimeType&& t1, corsika::units::si::TimeType&& t2)
+ : has_intersections_(true)
+ , intersections_(std::make_pair(t1, t2)) {}
+ Intersections(corsika::units::si::TimeType&& t)
+ : has_intersections_(true)
+ , intersections_(std::make_pair(
+ t,
+ std::numeric_limits<corsika::units::si::TimeType::value_type>::infinity() *
+ corsika::units::si::second)) {}
+
+ bool hasIntersections() const { return has_intersections_; }
+ ///! where did the trajectory currently enter the volume
+ corsika::units::si::TimeType getEntry() const { return intersections_.first; }
+ ///! where did the trajectory currently exit the volume
+ corsika::units::si::TimeType getExit() const { return intersections_.second; }
+
+ private:
+ bool has_intersections_;
+ std::pair<corsika::units::si::TimeType, corsika::units::si::TimeType> intersections_;
+ };
+
+} // namespace corsika::geometry
diff --git a/Framework/Geometry/Line.h b/Framework/Geometry/Line.h
index fd0835308..2d27e64df 100644
--- a/Framework/Geometry/Line.h
+++ b/Framework/Geometry/Line.h
@@ -14,19 +14,34 @@
namespace corsika::geometry {
+ /**
+ * \class Line
+ *
+ * A Line describes a movement in three dimensional space. It
+ * consists of a Point `$\vec{p_0}$` and and a speed-Vector
+ * `$\vec{v}$`, so that it can return GetPosition as
+ * `$\vec{p_0}*\vec{v}*t$` for any value of time `$t$`.
+ *
+ **/
+
class Line {
using VelocityVec = Vector<corsika::units::si::SpeedType::dimension_type>;
Point const r0;
VelocityVec const v0;
-
+
public:
+ Line() = delete;
+ Line(const Line&) = default;
+ Line(Line&&) = default;
+ Line& operator=(const Line&) = default;
Line(Point const& pR0, VelocityVec const& pV0)
: r0(pR0)
, v0(pV0) {}
Point GetPosition(corsika::units::si::TimeType t) const { return r0 + v0 * t; }
+ VelocityVec GetVelocity(corsika::units::si::TimeType) const { return v0; }
Point PositionFromArclength(corsika::units::si::LengthType l) const {
return r0 + v0.normalized() * l;
diff --git a/Framework/Geometry/Sphere.h b/Framework/Geometry/Sphere.h
index 52d3909c3..7d2ab39a3 100644
--- a/Framework/Geometry/Sphere.h
+++ b/Framework/Geometry/Sphere.h
@@ -10,6 +10,7 @@
#include <corsika/geometry/Point.h>
#include <corsika/geometry/Volume.h>
+#include <corsika/geometry/Line.h>
#include <corsika/units/PhysicalUnits.h>
namespace corsika::geometry {
@@ -28,8 +29,9 @@ namespace corsika::geometry {
return fRadius * fRadius > (fCenter - p).squaredNorm();
}
- auto& GetCenter() const { return fCenter; }
- auto GetRadius() const { return fRadius; }
+ const Point& GetCenter() const { return fCenter; }
+ LengthType GetRadius() const { return fRadius; }
+
};
} // namespace corsika::geometry
diff --git a/Framework/Geometry/Trajectory.h b/Framework/Geometry/Trajectory.h
index 1f9c24760..2c0e17afc 100644
--- a/Framework/Geometry/Trajectory.h
+++ b/Framework/Geometry/Trajectory.h
@@ -14,42 +14,192 @@
namespace corsika::geometry {
- template <typename T>
- class Trajectory : public T {
+ /**
+ * \class LineTrajectory
+ *
+ * A Trajectory is a description of a momvement of an object in
+ * three-dimensional space that describes the trajectory (connection
+ * between two Points in space), as well as the direction of motion
+ * at any given point.
+ *
+ * A Trajectory has a start `0` and an end `1`, where
+ * e.g. GetPosition(0) returns the start point and GetDirection(1)
+ * the direction of motion at the end. Values outside 0...1 are not
+ * defined.
+ *
+ * A Trajectory has a length in [m], GetLength, a duration in [s], GetDuration.
+ *
+ * Note: so far it is assumed that the speed (d|vec{r}|/dt) between
+ * start and end does not change and is constant for the entire
+ * Trajectory.
+ *
+ **/
- corsika::units::si::TimeType fTimeLength;
+ class LineTrajectory {
+
+ using VelocityVec = Vector<corsika::units::si::SpeedType::dimension_type>;
public:
- using T::ArcLength;
- using T::GetPosition;
+ LineTrajectory() = delete;
+ LineTrajectory(const LineTrajectory&) = default;
+ LineTrajectory(LineTrajectory&&) = default;
+ LineTrajectory& operator=(const LineTrajectory&) = default;
+ LineTrajectory(Line const& theLine, corsika::units::si::TimeType timeLength)
+ : line_(theLine)
+ , timeLength_(timeLength)
+ , timeStep_(timeLength)
+ , initialVelocity_(theLine.GetVelocity(timeLength * 0))
+ , finalVelocity_(theLine.GetVelocity(timeLength * 1)) {}
+ LineTrajectory(
+ Line const& theLine,
+ corsika::units::si::TimeType timeLength, // length of theLine (straight)
+ corsika::units::si::TimeType timeStep, // length of bend step (curved)
+ const VelocityVec& initialV, const VelocityVec& finalV)
+ : line_(theLine)
+ , timeLength_(timeLength)
+ , timeStep_(timeStep)
+ , initialVelocity_(initialV)
+ , finalVelocity_(finalV) {}
+
+ const Line& GetLine() const { return line_; }
+ Point GetPosition(double u) const { return line_.GetPosition(timeLength_ * u); }
+ VelocityVec GetVelocity(double u) const {
+ return initialVelocity_ * (1 - u) + finalVelocity_ * u;
+ }
+ Vector<corsika::units::si::dimensionless_d> GetDirection(double u) const {
+ return GetVelocity(u).normalized();
+ }
+
+ ///! duration along potentially bend trajectory
+ corsika::units::si::TimeType GetDuration(double u = 1) const { return u * timeStep_; }
+
+ ///! total length along potentially bend trajectory
+ corsika::units::si::LengthType GetLength(double u = 1) const {
+ using namespace corsika::units::si;
+ if (timeLength_ == 0_s) return 0_m;
+ return GetDistance(u) * timeStep_ / timeLength_;
+ }
+
+ ///! set new duration along potentially bend trajectory.
+ void SetLength(corsika::units::si::LengthType limit) {
+ SetDuration(line_.TimeFromArclength(limit));
+ }
+
+ ///! set new duration along potentially bend trajectory.
+ // Scale other properties by "limit/timeLength_"
+ void SetDuration(corsika::units::si::TimeType limit) {
+ using namespace corsika::units::si;
+ if (timeStep_ == 0_s) {
+ timeLength_ *= 0;
+ SetFinalVelocity(GetVelocity(0));
+ timeStep_ = limit;
+ } else {
+ const double scale = limit / timeStep_;
+ timeLength_ *= scale;
+ SetFinalVelocity(GetVelocity(scale));
+ timeStep_ = limit;
+ }
+ }
+
+ protected:
+ ///! total length along straight trajectory
+ corsika::units::si::LengthType GetDistance(double u) const {
+ assert(u <= 1);
+ assert(u >= 0);
+ return line_.ArcLength(0 * corsika::units::si::second, u * timeLength_);
+ }
+
+ void SetFinalVelocity(const VelocityVec& v) { finalVelocity_ = v; }
+
+ private:
+ Line line_;
+ corsika::units::si::TimeType timeLength_;
+ corsika::units::si::TimeType timeStep_;
+ VelocityVec initialVelocity_;
+ VelocityVec finalVelocity_;
+ };
+
+ /**
+ * \class LeapFrogTrajectory
+ *
+ *
+ **/
+
+ class LeapFrogTrajectory {
- Trajectory(T const& theT, corsika::units::si::TimeType timeLength)
- : T(theT)
- , fTimeLength(timeLength) {}
+ using VelocityVec = Vector<corsika::units::si::SpeedType::dimension_type>;
+ typedef corsika::geometry::Vector<corsika::units::si::magnetic_flux_density_d>
+ MagneticFieldVector;
- /*Point GetPosition(corsika::units::si::TimeType t) const {
- return fTraj.GetPosition(t + fTStart);
- }*/
+ public:
+ LeapFrogTrajectory() = delete;
+ LeapFrogTrajectory(const LeapFrogTrajectory&) = default;
+ LeapFrogTrajectory(LeapFrogTrajectory&&) = default;
+ LeapFrogTrajectory& operator=(const LeapFrogTrajectory&) = default;
+ LeapFrogTrajectory(const Point& pos, const VelocityVec& initialVelocity,
+ MagneticFieldVector Bfield,
+ const decltype(square(corsika::units::si::meter) /
+ (square(corsika::units::si::second) *
+ corsika::units::si::volt)) k,
+ corsika::units::si::TimeType timeStep) // leap-from total length
+ : initialPosition_(pos)
+ , initialVelocity_(initialVelocity)
+ , initialDirection_(initialVelocity.normalized())
+ , magneticfield_(Bfield)
+ , k_(k)
+ , timeStep_(timeStep) {}
- Point GetPosition(double u) const { return T::GetPosition(fTimeLength * u); }
+ const Line GetLine(double u) const { return Line(GetPosition(u), GetVelocity(u)); }
+ Point GetPosition(double u) const {
+ Point position = initialPosition_ + initialVelocity_ * timeStep_ * u / 2;
+ VelocityVec velocity =
+ initialVelocity_ + initialVelocity_.cross(magneticfield_) * timeStep_ * u * k_;
+ return position + velocity * timeStep_ * u / 2;
+ // auto steplength_true = steplength_true * (1.0 + double(direction.norm())) /
+ // 2;
+ }
+ VelocityVec GetVelocity(double u) const {
+ return initialVelocity_ +
+ initialVelocity_.cross(magneticfield_) * timeStep_ * u * k_;
+ }
+
+ Vector<corsika::units::si::dimensionless_d> GetDirection(double u) const {
+ return GetVelocity(u).normalized();
+ }
- corsika::units::si::TimeType GetDuration() const { return fTimeLength; }
- corsika::units::si::LengthType GetLength() const { return GetDistance(fTimeLength); }
+ ///! duration along potentially bend trajectory
+ corsika::units::si::TimeType GetDuration(double u = 1) const {
+ return u * timeStep_ * (double(GetVelocity(u).norm()/initialVelocity_.norm()) + 1.0) / 2;
+ }
- corsika::units::si::LengthType GetDistance(corsika::units::si::TimeType t) const {
- assert(t <= fTimeLength);
- assert(t >= 0 * corsika::units::si::second);
- return T::ArcLength(0 * corsika::units::si::second, t);
+ ///! total length along potentially bend trajectory
+ corsika::units::si::LengthType GetLength(double u = 1) const {
+ using namespace corsika::units::si;
+ return timeStep_ * initialVelocity_.norm() * u;
}
- void LimitEndTo(corsika::units::si::LengthType limit) {
- fTimeLength = T::TimeFromArclength(limit);
+ ///! set new duration along potentially bend trajectory.
+ void SetLength(corsika::units::si::LengthType limit) {
+ using namespace corsika::units::si;
+ if (initialVelocity_.norm() == 0_m / 1_s) SetDuration(0_s);
+ SetDuration(limit / initialVelocity_.norm());
}
- auto NormalizedDirection() const {
- static_assert(std::is_same_v<T, corsika::geometry::Line>);
- return T::GetV0().normalized();
+ ///! set new duration along potentially bend trajectory.
+ // Scale other properties by "limit/timeLength_"
+ void SetDuration(corsika::units::si::TimeType limit) {
+ using namespace corsika::units::si;
+ timeStep_ = limit;
}
+
+ private:
+ Point initialPosition_;
+ VelocityVec initialVelocity_;
+ geometry::Vector<corsika::units::si::dimensionless_d> initialDirection_;
+ MagneticFieldVector magneticfield_;
+ decltype(square(corsika::units::si::meter) /
+ (square(corsika::units::si::second) * corsika::units::si::volt)) k_;
+ corsika::units::si::TimeType timeStep_;
};
} // namespace corsika::geometry
diff --git a/Framework/Geometry/testGeometry.cc b/Framework/Geometry/testGeometry.cc
index 6211b3cc4..e522d0ab5 100644
--- a/Framework/Geometry/testGeometry.cc
+++ b/Framework/Geometry/testGeometry.cc
@@ -226,12 +226,10 @@ TEST_CASE("Trajectories") {
.magnitude() == Approx(0).margin(absMargin));
auto const t = 1_s;
- Trajectory<Line> base(line, t);
+ LineTrajectory base(line, t);
CHECK(line.GetPosition(t).GetCoordinates() == base.GetPosition(1.).GetCoordinates());
- CHECK(base.ArcLength(1_s, 2_s) / 1_m == Approx(3));
-
- CHECK((base.NormalizedDirection().GetComponents(rootCS) -
+ CHECK((base.GetVelocity(0).normalized().GetComponents(rootCS) -
QuantityVector<dimensionless_d>{1, 0, 0})
.eVector.norm() == Approx(0).margin(absMargin));
}
@@ -263,10 +261,5 @@ TEST_CASE("Trajectories") {
.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(0_s, 1_s) / 1_m == Approx(5));
}
}
diff --git a/Framework/Units/PhysicalUnits.h b/Framework/Units/PhysicalUnits.h
index 1f207abbc..36a835e7e 100644
--- a/Framework/Units/PhysicalUnits.h
+++ b/Framework/Units/PhysicalUnits.h
@@ -57,8 +57,6 @@ namespace corsika::units::si {
/// defining cross section as area
using sigma_d = phys::units::area_d;
- // constexpr quantity<area_d> barn{Rep(1e-28L) * square(meter)};
-
/// add the unit-types
using LengthType = phys::units::quantity<phys::units::length_d, double>;
using TimeType = phys::units::quantity<phys::units::time_interval_d, double>;
@@ -79,6 +77,8 @@ namespace corsika::units::si {
phys::units::quantity<phys::units::dimensions<0, 0, -1>, double>;
using InverseGrammageType =
phys::units::quantity<phys::units::dimensions<2, -1, 0>, double>;
+ using MagneticFluxType =
+ phys::units::quantity<phys::units::magnetic_flux_density_d, double>;
template <typename DimFrom, typename DimTo>
auto constexpr ConversionFactorHEPToSI() {
diff --git a/Processes/CMakeLists.txt b/Processes/CMakeLists.txt
index 4e9f59df3..63e66a331 100644
--- a/Processes/CMakeLists.txt
+++ b/Processes/CMakeLists.txt
@@ -3,7 +3,10 @@ add_subdirectory (AnalyticProcessors)
add_subdirectory (ExampleProcessors)
# tracking
+add_subdirectory (Tracking)
add_subdirectory (TrackingLine)
+add_subdirectory (TrackingLeapFrogStraight)
+add_subdirectory (TrackingLeapFrogCurved)
# hadron interaction models
add_subdirectory (Sibyll)
add_subdirectory (QGSJetII)
diff --git a/Processes/CONEXSourceCut/testCONEXSourceCut.cc b/Processes/CONEXSourceCut/testCONEXSourceCut.cc
index 09e905349..4c17d7fb4 100644
--- a/Processes/CONEXSourceCut/testCONEXSourceCut.cc
+++ b/Processes/CONEXSourceCut/testCONEXSourceCut.cc
@@ -12,7 +12,6 @@
#include <corsika/environment/LayeredSphericalAtmosphereBuilder.h>
#include <corsika/environment/MediumPropertyModel.h>
#include <corsika/environment/UniformMagneticField.h>
-#include <corsika/environment/UniformMediumType.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/RootCoordinateSystem.h>
@@ -63,31 +62,11 @@ TEST_CASE("CONEXSourceCut") {
{{particles::Code::Nitrogen, particles::Code::Oxygen},
{0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now
- builder.addExponentialLayer<
- environment::UniformMediumType<environment::UniformMagneticField<
- SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
- 1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km, environment::EMediumType::eAir,
- geometry::Vector(rootCS, 0_T, 0_T, 1_T));
- builder.addExponentialLayer<
- environment::UniformMediumType<environment::UniformMagneticField<
- SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
- 1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km, environment::EMediumType::eAir,
- geometry::Vector(rootCS, 0_T, 0_T, 1_T));
- builder.addExponentialLayer<
- environment::UniformMediumType<environment::UniformMagneticField<
- SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
- 1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km, environment::EMediumType::eAir,
- geometry::Vector(rootCS, 0_T, 0_T, 1_T));
- builder.addExponentialLayer<
- environment::UniformMediumType<environment::UniformMagneticField<
- SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
- 540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km, environment::EMediumType::eAir,
- geometry::Vector(rootCS, 0_T, 0_T, 1_T));
- builder.addLinearLayer<environment::UniformMediumType<
- environment::UniformMagneticField<HomogeneousMedium<setup::EnvironmentInterface>>>>(
- 1e9_cm, 112.8_km, environment::EMediumType::eAir,
- geometry::Vector(rootCS, 0_T, 0_T, 1_T));
-
+ builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km);
+ builder.addExponentialLayer(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km);
+ builder.addExponentialLayer(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km);
+ builder.addExponentialLayer(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km);
+ builder.addLinearLayer(1e9_cm, 112.8_km);
builder.assemble(env);
const HEPEnergyType E0 = 1_PeV;
@@ -112,7 +91,8 @@ TEST_CASE("CONEXSourceCut") {
[[maybe_unused]] process::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
corsika::process::conex_source_cut::CONEXSourceCut conex(
- center, showerAxis, t, injectionHeight, E0, particles::Code::Proton);
+ center, showerAxis, t, injectionHeight, E0,
+ particles::GetPDG(particles::Code::Proton));
HEPEnergyType const Eem{1_PeV};
auto const momentum = showerAxis.GetDirection() * Eem;
diff --git a/Processes/ObservationPlane/ObservationPlane.cc b/Processes/ObservationPlane/ObservationPlane.cc
index f624c8e3a..6153c6356 100644
--- a/Processes/ObservationPlane/ObservationPlane.cc
+++ b/Processes/ObservationPlane/ObservationPlane.cc
@@ -31,12 +31,13 @@ ObservationPlane::ObservationPlane(
corsika::process::EProcessReturn ObservationPlane::DoContinuous(
setup::Stack::ParticleType& particle, setup::Trajectory const& trajectory) {
TimeType const timeOfIntersection =
- (plane_.GetCenter() - trajectory.GetR0()).dot(plane_.GetNormal()) /
- trajectory.GetV0().dot(plane_.GetNormal());
+ (plane_.GetCenter() - trajectory.GetLine().GetR0()).dot(plane_.GetNormal()) /
+ trajectory.GetLine().GetV0().dot(plane_.GetNormal());
if (timeOfIntersection < TimeType::zero()) { return process::EProcessReturn::eOk; }
- if (plane_.IsAbove(trajectory.GetR0()) == plane_.IsAbove(trajectory.GetPosition(1))) {
+ if (plane_.IsAbove(trajectory.GetLine().GetR0()) ==
+ plane_.IsAbove(trajectory.GetPosition(1))) {
return process::EProcessReturn::eOk;
}
@@ -44,8 +45,8 @@ corsika::process::EProcessReturn ObservationPlane::DoContinuous(
auto const displacement = trajectory.GetPosition(1) - plane_.GetCenter();
outputStream_ << static_cast<int>(particles::GetPDG(particle.GetPID())) << ' '
- << energy / 1_eV << ' '
- << displacement.dot(xAxis_) / 1_m << ' ' << displacement.dot(yAxis_) / 1_m
+ << energy / 1_eV << ' ' << displacement.dot(xAxis_) / 1_m << ' '
+ << displacement.dot(yAxis_) / 1_m
<< (trajectory.GetPosition(1) - plane_.GetCenter()).norm() / 1_m
<< std::endl;
@@ -68,56 +69,69 @@ LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const& vPa
chargeNumber = corsika::particles::GetChargeNumber(vParticle.GetPID());
}
auto const* currentLogicalVolumeNode = vParticle.GetNode();
- auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(vParticle.GetPosition());
- auto direction = trajectory.GetV0().normalized();
-
- if (chargeNumber != 0 && abs(plane_.GetNormal().dot(trajectory.GetV0().cross(magneticfield))) * 1_s / 1_m / 1_T > 1e-6) {
+ auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(
+ vParticle.GetPosition());
+ auto direction = trajectory.GetLine().GetV0().normalized();
+
+ if (chargeNumber != 0 &&
+ abs(plane_.GetNormal().dot(trajectory.GetLine().GetV0().cross(magneticfield))) *
+ 1_s / 1_m / 1_T >
+ 1e-6) {
auto const* currentLogicalVolumeNode = vParticle.GetNode();
- auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(vParticle.GetPosition());
- auto k = chargeNumber * corsika::units::constants::c * 1_eV / (vParticle.GetMomentum().norm() * 1_V);
-
- if (direction.dot(plane_.GetNormal()) * direction.dot(plane_.GetNormal()) -
- (plane_.GetNormal().dot(trajectory.GetR0() - plane_.GetCenter()) *
- plane_.GetNormal().dot(direction.cross(magneticfield)) * 2 * k) < 0)
- {
- return std::numeric_limits<double>::infinity() * 1_m;
- }
-
- LengthType MaxStepLength1 =
- ( sqrt(direction.dot(plane_.GetNormal()) * direction.dot(plane_.GetNormal()) -
- (plane_.GetNormal().dot(trajectory.GetR0() - plane_.GetCenter()) *
- plane_.GetNormal().dot(direction.cross(magneticfield)) * 2 * k)) -
- direction.dot(plane_.GetNormal()) / direction.GetNorm() ) /
- (plane_.GetNormal().dot(direction.cross(magneticfield)) * k);
-
- LengthType MaxStepLength2 =
- ( - sqrt(direction.dot(plane_.GetNormal()) * direction.dot(plane_.GetNormal()) -
- (plane_.GetNormal().dot(trajectory.GetR0() - plane_.GetCenter()) *
- plane_.GetNormal().dot(direction.cross(magneticfield)) * 2 * k)) -
- direction.dot(plane_.GetNormal()) / direction.GetNorm() ) /
- (plane_.GetNormal().dot(direction.cross(magneticfield)) * k);
-
+ auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(
+ vParticle.GetPosition());
+ auto k = chargeNumber * corsika::units::constants::c * 1_eV /
+ (vParticle.GetMomentum().norm() * 1_V);
+
+ if (direction.dot(plane_.GetNormal()) * direction.dot(plane_.GetNormal()) -
+ (plane_.GetNormal().dot(trajectory.GetLine().GetR0() - plane_.GetCenter()) *
+ plane_.GetNormal().dot(direction.cross(magneticfield)) * 2 * k) <
+ 0) {
+ return std::numeric_limits<double>::infinity() * 1_m;
+ }
+
+ LengthType MaxStepLength1 =
+ (sqrt(direction.dot(plane_.GetNormal()) * direction.dot(plane_.GetNormal()) -
+ (plane_.GetNormal().dot(trajectory.GetLine().GetR0() - plane_.GetCenter()) *
+ plane_.GetNormal().dot(direction.cross(magneticfield)) * 2 * k)) -
+ direction.dot(plane_.GetNormal()) / direction.GetNorm()) /
+ (plane_.GetNormal().dot(direction.cross(magneticfield)) * k);
+
+ LengthType MaxStepLength2 =
+ (-sqrt(
+ direction.dot(plane_.GetNormal()) * direction.dot(plane_.GetNormal()) -
+ (plane_.GetNormal().dot(trajectory.GetLine().GetR0() - plane_.GetCenter()) *
+ plane_.GetNormal().dot(direction.cross(magneticfield)) * 2 * k)) -
+ direction.dot(plane_.GetNormal()) / direction.GetNorm()) /
+ (plane_.GetNormal().dot(direction.cross(magneticfield)) * k);
+
if (MaxStepLength1 <= 0_m && MaxStepLength2 <= 0_m) {
return std::numeric_limits<double>::infinity() * 1_m;
} else if (MaxStepLength1 <= 0_m || MaxStepLength2 < MaxStepLength1) {
- std::cout << " steplength to obs plane 2: " << MaxStepLength2 << std::endl;
- return MaxStepLength2 * (direction + direction.cross(magneticfield) * MaxStepLength2 * k / 2).norm() * 1.001;
+ std::cout << " steplength to obs plane 2: " << MaxStepLength2 << std::endl;
+ return MaxStepLength2 *
+ (direction + direction.cross(magneticfield) * MaxStepLength2 * k / 2)
+ .norm() *
+ 1.001;
} else if (MaxStepLength2 <= 0_m || MaxStepLength1 < MaxStepLength2) {
- std::cout << " steplength to obs plane 1: " << MaxStepLength1 << std::endl;
- return MaxStepLength1 * (direction + direction.cross(magneticfield) * MaxStepLength2 * k / 2).norm() * 1.001;
+ std::cout << " steplength to obs plane 1: " << MaxStepLength1 << std::endl;
+ return MaxStepLength1 *
+ (direction + direction.cross(magneticfield) * MaxStepLength2 * k / 2)
+ .norm() *
+ 1.001;
}
- }
+ }
TimeType const timeOfIntersection =
- (plane_.GetCenter() - trajectory.GetR0()).dot(plane_.GetNormal()) /
- trajectory.GetV0().dot(plane_.GetNormal());
+ (plane_.GetCenter() - trajectory.GetLine().GetR0()).dot(plane_.GetNormal()) /
+ trajectory.GetLine().GetV0().dot(plane_.GetNormal());
if (timeOfIntersection < TimeType::zero()) {
return std::numeric_limits<double>::infinity() * 1_m;
}
- auto const pointOfIntersection = trajectory.GetPosition(timeOfIntersection);
+ auto const pointOfIntersection = trajectory.GetLine().GetPosition(timeOfIntersection);
std::cout << " obs plane non b-field " << std::endl;
- return (trajectory.GetR0() - pointOfIntersection).norm() * 1.0001;
+ return (trajectory.GetLine().GetR0() - pointOfIntersection).norm() * 1.0001;
}
void ObservationPlane::ShowResults() const {
diff --git a/Processes/ObservationPlane/testObservationPlane.cc b/Processes/ObservationPlane/testObservationPlane.cc
index 42837e8b2..3010d34b1 100644
--- a/Processes/ObservationPlane/testObservationPlane.cc
+++ b/Processes/ObservationPlane/testObservationPlane.cc
@@ -51,7 +51,7 @@ TEST_CASE("ContinuousProcess interface", "[proccesses][observation_plane]") {
Vector<units::si::SpeedType::dimension_type> vec(cs, 0_m / second, 0_m / second,
-units::constants::c);
Line line(start, vec);
- Trajectory<Line> track(line, 12_m / units::constants::c);
+ LineTrajectory track(line, 12_m / units::constants::c);
particle.SetPosition(Point(cs, {1_m, 1_m, 10_m})); // moving already along -z
diff --git a/Processes/StackInspector/testStackInspector.cc b/Processes/StackInspector/testStackInspector.cc
index 9814efa02..3d4cdec2c 100644
--- a/Processes/StackInspector/testStackInspector.cc
+++ b/Processes/StackInspector/testStackInspector.cc
@@ -31,7 +31,7 @@ TEST_CASE("StackInspector", "[processes]") {
geometry::Vector<units::si::SpeedType::dimension_type> v(rootCS, 0_m / second,
0_m / second, 1_m / second);
geometry::Line line(origin, v);
- geometry::Trajectory<geometry::Line> track(line, 10_s);
+ geometry::LineTrajectory track(line, 10_s);
TestCascadeStack stack;
stack.Clear();
diff --git a/Processes/Tracking/CMakeLists.txt b/Processes/Tracking/CMakeLists.txt
new file mode 100644
index 000000000..d0aebb6d3
--- /dev/null
+++ b/Processes/Tracking/CMakeLists.txt
@@ -0,0 +1,34 @@
+set (
+ MODEL_HEADERS
+ Intersect.hpp
+ )
+
+set (
+ MODEL_NAMESPACE
+ corsika/process/tracking
+ )
+
+add_library (ProcessTrackingIntersects INTERFACE)
+CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessTrackingIntersects ${MODEL_NAMESPACE} ${MODEL_HEADERS})
+
+# target dependencies on other libraries (also the header onlys)
+target_link_libraries (
+ ProcessTrackingIntersects
+ INTERFACE
+ CORSIKAsetup
+ CORSIKAutilities
+ CORSIKAenvironment
+ CORSIKAunits
+ CORSIKAgeometry
+ CORSIKAlogging
+ )
+
+target_include_directories (
+ ProcessTrackingIntersects
+ INTERFACE
+ $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
+ $<INSTALL_INTERFACE:include/include>
+ )
+
+install (FILES ${MODEL_HEADERS} DESTINATION include/${MODEL_NAMESPACE})
+
diff --git a/Processes/Tracking/Intersect.hpp b/Processes/Tracking/Intersect.hpp
new file mode 100644
index 000000000..19e4819bf
--- /dev/null
+++ b/Processes/Tracking/Intersect.hpp
@@ -0,0 +1,128 @@
+#pragma once
+
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/Vector.h>
+#include <corsika/units/PhysicalUnits.h>
+#include <corsika/logging/Logging.h>
+#include <corsika/geometry/Intersections.hpp>
+
+namespace corsika::process::tracking {
+
+ /**
+ * \class Intersect
+ *
+ *
+ *
+ **/
+
+ template <typename TDerived>
+ class Intersect {
+
+ protected:
+ template <typename TParticle>
+ auto nextIntersect(const TParticle& particle) const {
+ using namespace corsika::units::si;
+ using namespace corsika::geometry;
+
+ const Point& initialPosition = particle.GetPosition();
+
+ typedef
+ typename std::remove_reference<decltype(*particle.GetNode())>::type node_type;
+ node_type& volumeNode = *particle.GetNode();
+ C8LOG_DEBUG("volumeNode={}, numericallyInside={} ", fmt::ptr(&volumeNode),
+ volumeNode.GetVolume().Contains(initialPosition));
+
+ auto const velocity =
+ particle.GetMomentum() / particle.GetEnergy() * corsika::units::constants::c;
+
+ // for the event of magnetic fields and curved trajectories, we need to limit
+ // maximum step-length since we need to follow curved
+ // trajectories segment-wise -- at least if we don't employ concepts as "Helix
+ // Trajectories" or similar
+ const auto& magneticfield =
+ volumeNode.GetModelProperties().GetMagneticField(initialPosition);
+ const auto magnitudeB = magneticfield.norm();
+ const int chargeNumber = particle.GetChargeNumber();
+ auto const momentumVerticalMag =
+ particle.GetMomentum() -
+ particle.GetMomentum().parallelProjectionOnto(magneticfield);
+ LengthType const gyroradius =
+ (chargeNumber == 0 || magnitudeB == 0_T
+ ? std::numeric_limits<TimeType::value_type>::infinity() * 1_m
+ : momentumVerticalMag.norm() * 1_V /
+ (corsika::units::constants::c * abs(chargeNumber) * magnitudeB *
+ 1_eV));
+ const double maxRadians = 0.01;
+ const LengthType steplimit = 2 * cos(maxRadians) * sin(maxRadians) * gyroradius;
+ C8LOG_DEBUG("gyroradius {}, steplimit: {} m = {} s", gyroradius, steplimit,
+ steplimit / velocity.norm());
+
+ // start values:
+ TimeType minTime = steplimit / velocity.norm();
+ node_type* minNode = &volumeNode;
+
+ // determine the first geometric collision with any other Volume boundary
+
+ // first check, where we leave the current volume
+ // this assumes our convention, that all Volume primitives must be convex
+ // thus, the last entry is always the exit point
+ const Intersections time_intersections_curr =
+ TDerived::Intersect(particle, volumeNode);
+ C8LOG_TRACE("curr node {}, parent node {} ", fmt::ptr(&volumeNode),
+ fmt::ptr(volumeNode.GetParent()));
+ C8LOG_DEBUG("intersection times with currentLogicalVolumeNode: {} s and {} s",
+ time_intersections_curr.getEntry() / 1_s,
+ time_intersections_curr.getExit() / 1_s);
+ if (time_intersections_curr.getExit() <= minTime) {
+ minTime =
+ time_intersections_curr.getExit(); // we exit currentLogicalVolumeNode here
+ minNode = volumeNode.GetParent();
+ }
+
+ // where do we collide with any of the next-tree-level volumes
+ // entirely contained by currentLogicalVolumeNode
+ for (const auto& node : volumeNode.GetChildNodes()) {
+
+ const Intersections time_intersections = TDerived::Intersect(particle, *node);
+ C8LOG_DEBUG("intersection times with child volume {} : enter {} s, exit {} s",
+ fmt::ptr(node), time_intersections.getEntry() / 1_s,
+ time_intersections.getExit() / 1_s);
+
+ const auto t_entry = time_intersections.getEntry();
+ const auto t_exit = time_intersections.getExit();
+ C8LOG_TRACE("children t-entry: {}, t-exit: {}, smaller? {} ", t_entry, t_exit,
+ t_entry <= minTime);
+ // note, theoretically t can even be smaller than 0 since we
+ // KNOW we can't yet be in this volume yet, so we HAVE TO
+ // enter it IF exit point is not also in the "past"!
+ if (t_exit > 0_s && t_entry <= minTime) { // enter volumen child here
+ minTime = t_entry;
+ minNode = node.get();
+ }
+ }
+
+ // these are volumes from the previous tree-level that are cut-out partly from the
+ // current volume
+ for (node_type* node : volumeNode.GetExcludedNodes()) {
+
+ const Intersections time_intersections = TDerived::Intersect(particle, *node);
+ C8LOG_DEBUG("intersection times with exclusion volume {} : enter {} s, exit {} s",
+ fmt::ptr(node), time_intersections.getEntry() / 1_s,
+ time_intersections.getExit() / 1_s);
+ const auto t_entry = time_intersections.getEntry();
+ const auto t_exit = time_intersections.getExit();
+ C8LOG_TRACE("children t-entry: {}, t-exit: {}, smaller? {} ", t_entry, t_exit,
+ t_entry <= minTime);
+ // note, theoretically t can even be smaller than 0 since we
+ // KNOW we can't yet be in this volume yet, so we HAVE TO
+ // enter it IF exit point is not also in the "past"!
+ if (t_exit > 0_s && t_entry <= minTime) { // enter volumen child here
+ minTime = t_entry;
+ minNode = node;
+ }
+ }
+ C8LOG_TRACE("t-intersect: {}, node {} ", minTime, fmt::ptr(minNode));
+ return std::make_tuple(minTime, minNode);
+ }
+ }; // namespace corsika::process::tracking
+} // namespace corsika::process::tracking
diff --git a/Processes/TrackingLeapFrogCurved/CMakeLists.txt b/Processes/TrackingLeapFrogCurved/CMakeLists.txt
new file mode 100644
index 000000000..141929f27
--- /dev/null
+++ b/Processes/TrackingLeapFrogCurved/CMakeLists.txt
@@ -0,0 +1,44 @@
+set (
+ MODEL_HEADERS
+ Tracking.h
+ )
+
+set (
+ MODEL_NAMESPACE
+ corsika/process/tracking_leapfrog_curved
+ )
+
+add_library (ProcessTrackingLeapFrogCurved INTERFACE)
+CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessTrackingLeapFrogCurved ${MODEL_NAMESPACE} ${MODEL_HEADERS})
+
+# target dependencies on other libraries (also the header onlys)
+target_link_libraries (
+ ProcessTrackingLeapFrogCurved
+ INTERFACE
+ ProcessTrackingIntersects
+ CORSIKAsetup
+ CORSIKAutilities
+ CORSIKAenvironment
+ CORSIKAunits
+ CORSIKAenvironment
+ CORSIKAgeometry
+ CORSIKAlogging
+ )
+
+target_include_directories (
+ ProcessTrackingLeapFrogCurved
+ INTERFACE
+ $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
+ $<INSTALL_INTERFACE:include/include>
+ )
+
+install (FILES ${MODEL_HEADERS} DESTINATION include/${MODEL_NAMESPACE})
+
+# #-- -- -- -- -- -- -- -- -- --
+# #code unit testing
+CORSIKA_ADD_TEST (testTrackingLeapFrogCurved)
+target_link_libraries (
+ testTrackingLeapFrogCurved
+ ProcessTrackingLeapFrogCurved
+ CORSIKAtesting
+)
diff --git a/Processes/TrackingLeapFrogCurved/Tracking.h b/Processes/TrackingLeapFrogCurved/Tracking.h
new file mode 100644
index 000000000..f97ec44e0
--- /dev/null
+++ b/Processes/TrackingLeapFrogCurved/Tracking.h
@@ -0,0 +1,205 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#pragma once
+
+#include <corsika/process/tracking_line/Tracking.h>
+#include <corsika/process/tracking/Intersect.hpp>
+#include <corsika/geometry/Line.h>
+#include <corsika/geometry/Plane.h>
+#include <corsika/geometry/Sphere.h>
+#include <corsika/geometry/Trajectory.h>
+#include <corsika/geometry/Vector.h>
+#include <corsika/geometry/Intersections.hpp>
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/units/PhysicalUnits.h>
+#include <corsika/utl/quartic.h>
+#include <corsika/logging/Logging.h>
+
+#include <type_traits>
+#include <utility>
+
+#include <fstream>
+
+namespace corsika::process {
+
+ namespace tracking_leapfrog_curved {
+
+ typedef corsika::geometry::Vector<corsika::units::si::magnetic_flux_density_d>
+ MagneticFieldVector;
+
+ /**
+ * \function LeapFrogStep
+ *
+ * Performs one leap-frog step consistent of two halve-steps with steplength/2
+ **/
+ template <typename TParticle>
+ auto LeapFrogStep(const TParticle& particle,
+ corsika::units::si::LengthType steplength) {
+ using namespace corsika::units::si;
+ if (particle.GetMomentum().norm() == 0_GeV) {
+ return std::make_tuple(particle.GetPosition(), particle.GetMomentum() / 1_GeV,
+ double(0));
+ } // charge of the particle
+ const int chargeNumber = particle.GetChargeNumber();
+ auto const* currentLogicalVolumeNode = particle.GetNode();
+ auto magneticfield =
+ currentLogicalVolumeNode->GetModelProperties().GetMagneticField(
+ particle.GetPosition());
+ geometry::Vector<SpeedType::dimension_type> velocity =
+ particle.GetMomentum() / particle.GetEnergy() * corsika::units::constants::c;
+ decltype(corsika::units::si::meter /
+ (corsika::units::si::second * corsika::units::si::volt)) k =
+ chargeNumber * corsika::units::constants::cSquared * 1_eV /
+ (velocity.norm() * particle.GetEnergy() * 1_V);
+ geometry::Vector<dimensionless_d> direction = velocity.normalized();
+ auto position = particle.GetPosition(); // First Movement
+ // assuming magnetic field does not change during movement
+ position =
+ position + direction * steplength / 2; // Change of direction by magnetic field
+ direction =
+ direction + direction.cross(magneticfield) * steplength * k; // Second Movement
+ position = position + direction * steplength / 2;
+ auto steplength_true = steplength_true * (1.0 + double(direction.norm())) / 2;
+ return std::make_tuple(position, direction.normalized(), steplength_true);
+ }
+
+ /**
+ * \class Tracking
+ *
+ * The class tracking_leapfrog_curved::Tracking is based on the
+ * Bachelor thesis of Andre Schmidt (KIT). It implements a
+ * two-step leap-frog algorithm, but with analytically exact geometric
+ * intersections between leap-frog steps and geometric volumes
+ * (spheres, planes).
+ *
+ **/
+
+ class Tracking : public corsika::process::tracking::Intersect<Tracking> {
+
+ public:
+ template <typename TParticle>
+ auto GetTrack(TParticle const& particle) {
+ using namespace corsika::units::si;
+ using namespace corsika::geometry;
+ geometry::Vector<SpeedType::dimension_type> const initialVelocity =
+ particle.GetMomentum() / particle.GetEnergy() * corsika::units::constants::c;
+
+ auto const position = particle.GetPosition();
+ C8LOG_DEBUG(
+ "Tracking pid: {}"
+ " , E = {} GeV",
+ particle.GetPID(), particle.GetEnergy() / 1_GeV);
+ C8LOG_DEBUG("Tracking pos: {}", position.GetCoordinates());
+ C8LOG_DEBUG("Tracking E: {} GeV", particle.GetEnergy() / 1_GeV);
+ C8LOG_DEBUG("Tracking p: {} GeV",
+ particle.GetMomentum().GetComponents() / 1_GeV);
+ C8LOG_DEBUG("Tracking v: {} ", initialVelocity.GetComponents());
+
+ typedef
+ typename std::remove_reference<decltype(*particle.GetNode())>::type node_type;
+ node_type& volumeNode = *particle.GetNode();
+
+ // traverse the environment volume tree and find next
+ // intersection
+ auto [minTime, minNode] = tracking::Intersect<Tracking>::nextIntersect(particle);
+
+ const int chargeNumber = particle.GetChargeNumber();
+ const auto& magneticfield =
+ volumeNode.GetModelProperties().GetMagneticField(position);
+ const auto k = chargeNumber * corsika::units::constants::cSquared * 1_eV /
+ (particle.GetEnergy() * 1_V);
+ return std::make_tuple(
+ geometry::LeapFrogTrajectory(position, initialVelocity, magneticfield, k,
+ minTime), // trajectory
+ minNode); // next volume node
+ }
+
+ template <typename TParticle, typename TMedium>
+ static geometry::Intersections Intersect(const TParticle& particle,
+ const corsika::geometry::Sphere& sphere,
+ const TMedium& medium) {
+ using namespace corsika::units::si;
+ const int chargeNumber = particle.GetChargeNumber();
+ const auto& position = particle.GetPosition();
+ const auto& magneticfield = medium.GetMagneticField(position);
+
+ if (chargeNumber == 0 || magneticfield.norm() == 0_T) {
+ return tracking_line::Tracking::Intersect(particle, sphere, medium);
+ }
+
+ const geometry::Vector<SpeedType::dimension_type> velocity =
+ particle.GetMomentum() / particle.GetEnergy() * corsika::units::constants::c;
+ const auto absVelocity = velocity.norm();
+ auto energy = particle.GetEnergy();
+ auto k = chargeNumber * corsika::units::constants::cSquared * 1_eV /
+ (absVelocity * energy * 1_V);
+ const geometry::Vector<dimensionless_d> directionBefore =
+ velocity.normalized(); // determine steplength to next volume
+
+ const double a =
+ ((directionBefore.cross(magneticfield)).dot(position - sphere.GetCenter()) *
+ k +
+ 1) *
+ 4 /
+ (1_m * 1_m * (directionBefore.cross(magneticfield)).GetSquaredNorm() * k * k);
+ const double b = directionBefore.dot(position - sphere.GetCenter()) * 8 /
+ ((directionBefore.cross(magneticfield)).GetSquaredNorm() * k *
+ k * 1_m * 1_m * 1_m);
+ const double c = ((position - sphere.GetCenter()).GetSquaredNorm() -
+ (sphere.GetRadius() * sphere.GetRadius())) *
+ 4 /
+ ((directionBefore.cross(magneticfield)).GetSquaredNorm() * k *
+ k * 1_m * 1_m * 1_m * 1_m);
+ std::complex<double>* solutions = solve_quartic(0, a, b, c);
+ LengthType d_enter, d_exit;
+ int first = 0;
+ for (int i = 0; i < 4; i++) {
+ if (solutions[i].imag() == 0) {
+ LengthType time = solutions[i].real() * 1_m;
+ C8LOG_TRACE("Solutions for current Volume: {} ", time);
+ if (first == 0) {
+ d_enter = time;
+ } else {
+ if (time < d_enter) {
+ d_exit = d_enter;
+ d_enter = time;
+ } else {
+ d_exit = time;
+ }
+ }
+ first++;
+ }
+ }
+ delete[] solutions;
+
+ if (first != 2) {
+ C8LOG_DEBUG("no intersection! count={}", first);
+ return geometry::Intersections();
+ }
+ return geometry::Intersections(d_enter / absVelocity, d_exit / absVelocity);
+ }
+
+ template <typename TParticle, typename TBaseNodeType>
+ static geometry::Intersections Intersect(const TParticle& particle,
+ const TBaseNodeType& volumeNode) {
+ const geometry::Sphere* sphere =
+ dynamic_cast<const geometry::Sphere*>(&volumeNode.GetVolume());
+ if (sphere) {
+ return Intersect(particle, *sphere, volumeNode.GetModelProperties());
+ }
+ throw std::runtime_error(
+ "The Volume type provided is not supported in Intersect(particle, node)");
+ }
+ };
+
+ } // namespace tracking_leapfrog_curved
+
+} // namespace corsika::process
diff --git a/Processes/TrackingLeapFrogCurved/testTrackingLeapFrogCurved.cc b/Processes/TrackingLeapFrogCurved/testTrackingLeapFrogCurved.cc
new file mode 100644
index 000000000..cf45ccc48
--- /dev/null
+++ b/Processes/TrackingLeapFrogCurved/testTrackingLeapFrogCurved.cc
@@ -0,0 +1,127 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#include <corsika/process/tracking_leapfrog_curved/Tracking.h>
+
+#include <corsika/particles/ParticleProperties.h>
+
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/Sphere.h>
+#include <corsika/geometry/Vector.h>
+
+#include <corsika/setup/SetupEnvironment.h>
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+using corsika::setup::Trajectory;
+
+#include <catch2/catch.hpp>
+
+using namespace corsika;
+using namespace corsika::particles;
+using namespace corsika::process;
+using namespace corsika::units;
+using namespace corsika::geometry;
+using namespace corsika::units::si;
+
+template <typename T>
+int sgn(T val) {
+ return (T(0) < val) - (val < T(0));
+}
+
+TEST_CASE("TrackingLeapfrog_Curved") {
+
+ logging::SetLevel(logging::level::trace);
+
+ const HEPEnergyType P0 = 10_GeV;
+
+ auto PID = GENERATE(as<Code>{}, Code::MuPlus, Code::MuPlus, Code::Gamma);
+ auto Bfield = GENERATE(as<MagneticFluxType>{}, 0_T, 50_uT, -50_uT);
+ auto outer = GENERATE(as<bool>{}, true, false);
+
+ SECTION(fmt::format("Tracking PID={}, Bfield={} uT, from outside={}", PID,
+ Bfield / 1_uT, outer)) {
+
+ C8LOG_DEBUG(
+ "********************\n TEST section PID={}, Bfield={} "
+ "uT, outer (?)={}",
+ PID, Bfield / 1_uT, outer);
+
+ const int chargeNumber = GetChargeNumber(PID);
+ LengthType radius = 10_m;
+ int deflect = 0;
+ if (chargeNumber != 0 and Bfield != 0_T) {
+ deflect = -sgn(chargeNumber) * sgn(Bfield / 1_T); // direction of deflection
+ LengthType const gyroradius =
+ P0 * 1_V / (constants::c * abs(chargeNumber) * abs(Bfield) * 1_eV);
+ radius = gyroradius;
+ }
+
+ auto [env, csPtr, worldPtr] = setup::testing::setupEnvironment(Code::Oxygen, Bfield);
+ { [[maybe_unused]] const auto& env_dummy = env; }
+ auto const& cs = *csPtr;
+
+ tracking_leapfrog_curved::Tracking tracking;
+ using tracking_leapfrog_curved::MagneticFieldVector;
+ Point const center(cs, {0_m, 0_m, 0_m});
+ auto target = setup::Environment::CreateNode<geometry::Sphere>(center, radius);
+
+ using MyHomogeneousModel =
+ environment::MediumPropertyModel<environment::UniformMagneticField<
+ environment::HomogeneousMedium<setup::EnvironmentInterface>>>;
+
+ MagneticFieldVector magneticfield(cs, 0_T, 0_T, Bfield);
+ target->SetModelProperties<MyHomogeneousModel>(
+ environment::Medium::AirDry1Atm, magneticfield, 1_g / (1_m * 1_m * 1_m),
+ environment::NuclearComposition(std::vector<Code>{Code::Oxygen},
+ std::vector<float>{1.}));
+ auto* targetPtr = target.get();
+ worldPtr->AddChild(std::move(target));
+
+ auto [stack, viewPtr] = setup::testing::setupStack(PID, 0, 0, P0, targetPtr, cs);
+ { [[maybe_unused]] auto& viewPtr_dum = viewPtr; }
+ auto particle = stack->first();
+ // Note: momentum in X-direction
+ // magnetic field in Z-direction
+ // put particle on x_start, 0, 0
+ // expect intersections somewere in +-y_start
+
+ if (outer) {
+ particle.SetNode(worldPtr); // set particle inside "target" volume
+ } else {
+ particle.SetNode(targetPtr); // set particle outside "target" volume
+ }
+ particle.SetPosition(Point(cs, -radius, 0_m, 0_m));
+
+ auto [traj, nextVol] = tracking.GetTrack(particle);
+ particle.SetNode(nextVol);
+ particle.SetPosition(traj.GetPosition(1));
+ particle.SetMomentum(traj.GetDirection(1) * particle.GetMomentum().norm());
+ if (outer) {
+ // now we know we are in target volume, depending on "outer"
+ CHECK(traj.GetLength(1) == 0_m);
+ CHECK(nextVol == targetPtr);
+ }
+ // move forward, until we leave target volume
+ while (nextVol == targetPtr) {
+ const auto [traj2, nextVol2] = tracking.GetTrack(particle);
+ nextVol = nextVol2;
+ particle.SetNode(nextVol);
+ particle.SetPosition(traj2.GetPosition(1));
+ particle.SetMomentum(traj2.GetDirection(1) * particle.GetMomentum().norm());
+ }
+ CHECK(nextVol == worldPtr);
+
+ Point pointCheck(cs, (deflect == 0 ? radius : 0_m), (deflect * radius), 0_m);
+
+ C8LOG_DEBUG("testTrackingLineStack: deflect={}, momentum={}, pos={}, pos_check={}",
+ deflect, particle.GetMomentum().GetComponents(),
+ particle.GetPosition().GetCoordinates(), pointCheck.GetCoordinates());
+
+ CHECK((particle.GetPosition() - pointCheck).norm() / radius == Approx(0).margin(1e-3));
+ }
+}
diff --git a/Processes/TrackingLeapFrogCurved/testTrackingLeapFrogStraight.cc b/Processes/TrackingLeapFrogCurved/testTrackingLeapFrogStraight.cc
new file mode 100644
index 000000000..c96c2ffc8
--- /dev/null
+++ b/Processes/TrackingLeapFrogCurved/testTrackingLeapFrogStraight.cc
@@ -0,0 +1,127 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#include <corsika/process/tracking_bfield/Tracking.h>
+
+#include <corsika/particles/ParticleProperties.h>
+
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/Sphere.h>
+#include <corsika/geometry/Vector.h>
+
+#include <corsika/setup/SetupEnvironment.h>
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+using corsika::setup::Trajectory;
+
+#include <catch2/catch.hpp>
+
+using namespace corsika;
+using namespace corsika::particles;
+using namespace corsika::process;
+using namespace corsika::units;
+using namespace corsika::geometry;
+using namespace corsika::units::si;
+
+template <typename T>
+int sgn(T val) {
+ return (T(0) < val) - (val < T(0));
+}
+
+TEST_CASE("TrackingBField") {
+
+ logging::SetLevel(logging::level::trace);
+
+ const HEPEnergyType P0 = 10_GeV;
+
+ auto PID = GENERATE(as<Code>{}, Code::MuPlus, Code::MuPlus, Code::Gamma);
+ auto Bfield = GENERATE(as<MagneticFluxType>{}, 0_T, 50_uT, -50_uT);
+ auto outer = GENERATE(as<bool>{}, true, false);
+
+ SECTION(fmt::format("Tracking PID={}, Bfield={} uT, from outside={}", PID,
+ Bfield / 1_uT, outer)) {
+
+ C8LOG_DEBUG(
+ "********************\n TEST section PID={}, Bfield={} "
+ "uT, outer (?)={}",
+ PID, Bfield / 1_uT, outer);
+
+ const int chargeNumber = GetChargeNumber(PID);
+ LengthType radius = 10_m;
+ int deflect = 0;
+ if (chargeNumber != 0 and Bfield != 0_T) {
+ deflect = -sgn(chargeNumber) * sgn(Bfield / 1_T); // direction of deflection
+ LengthType const gyroradius =
+ P0 * 1_V / (constants::c * abs(chargeNumber) * abs(Bfield) * 1_eV);
+ radius = gyroradius;
+ }
+
+ auto [env, csPtr, worldPtr] = setup::testing::setupEnvironment(Code::Oxygen, Bfield);
+ { [[maybe_unused]] const auto& env_dummy = env; }
+ auto const& cs = *csPtr;
+
+ tracking_bfield::Tracking tracking;
+ using tracking_bfield::MagneticFieldVector;
+ Point const center(cs, {0_m, 0_m, 0_m});
+ auto target = setup::Environment::CreateNode<geometry::Sphere>(center, radius);
+
+ using MyHomogeneousModel =
+ environment::MediumPropertyModel<environment::UniformMagneticField<
+ environment::HomogeneousMedium<setup::EnvironmentInterface>>>;
+
+ MagneticFieldVector magneticfield(cs, 0_T, 0_T, Bfield);
+ target->SetModelProperties<MyHomogeneousModel>(
+ environment::Medium::AirDry1Atm, magneticfield, 1_g / (1_m * 1_m * 1_m),
+ environment::NuclearComposition(std::vector<Code>{Code::Oxygen},
+ std::vector<float>{1.}));
+ auto* targetPtr = target.get();
+ worldPtr->AddChild(std::move(target));
+
+ auto [stack, viewPtr] = setup::testing::setupStack(PID, 0, 0, P0, targetPtr, cs);
+ { [[maybe_unused]] auto& viewPtr_dum = viewPtr; }
+ auto particle = stack->first();
+ // Note: momentum in X-direction
+ // magnetic field in Z-direction
+ // put particle on x_start, 0, 0
+ // expect intersections somewere in +-y_start
+
+ if (outer) {
+ particle.SetNode(worldPtr); // set particle inside "target" volume
+ } else {
+ particle.SetNode(targetPtr); // set particle outside "target" volume
+ }
+ particle.SetPosition(Point(cs, -radius, 0_m, 0_m));
+
+ auto [traj, nextVol] = tracking.GetTrack(particle);
+ particle.SetNode(nextVol);
+ particle.SetPosition(traj.GetPosition(1));
+ particle.SetMomentum(traj.GetDirection(1) * particle.GetMomentum().norm());
+ if (outer) {
+ // now we know we are in target volume, depending on "outer"
+ CHECK(traj.GetLength(1) == 0_m);
+ CHECK(nextVol == targetPtr);
+ }
+ // move forward, until we leave target volume
+ while (nextVol == targetPtr) {
+ const auto [traj2, nextVol2] = tracking.GetTrack(particle);
+ nextVol = nextVol2;
+ particle.SetNode(nextVol);
+ particle.SetPosition(traj2.GetPosition(1));
+ particle.SetMomentum(traj2.GetDirection(1) * particle.GetMomentum().norm());
+ }
+ CHECK(nextVol == worldPtr);
+
+ Point pointCheck(cs, (deflect == 0 ? radius : 0_m), (deflect * radius), 0_m);
+
+ C8LOG_DEBUG("testTrackingLineStack: deflect={}, momentum={}, pos={}, pos_check={}",
+ deflect, particle.GetMomentum().GetComponents(),
+ particle.GetPosition().GetCoordinates(), pointCheck.GetCoordinates());
+
+ CHECK((particle.GetPosition() - pointCheck).norm() / radius == Approx(0).margin(1e-3));
+ }
+}
diff --git a/Processes/TrackingLeapFrogStraight/CMakeLists.txt b/Processes/TrackingLeapFrogStraight/CMakeLists.txt
new file mode 100644
index 000000000..afd5d99c7
--- /dev/null
+++ b/Processes/TrackingLeapFrogStraight/CMakeLists.txt
@@ -0,0 +1,43 @@
+set (
+ MODEL_HEADERS
+ Tracking.h
+ )
+
+set (
+ MODEL_NAMESPACE
+ corsika/process/tracking_leapfrag_straight
+ )
+
+add_library (ProcessTrackingLeapFrogStraight INTERFACE)
+CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessTrackingLeapFrogStraight ${MODEL_NAMESPACE} ${MODEL_HEADERS})
+
+# target dependencies on other libraries (also the header onlys)
+target_link_libraries (
+ ProcessTrackingLeapFrogStraight
+ INTERFACE
+ CORSIKAsetup
+ CORSIKAutilities
+ CORSIKAenvironment
+ CORSIKAunits
+ CORSIKAenvironment
+ CORSIKAgeometry
+ CORSIKAlogging
+ )
+
+target_include_directories (
+ ProcessTrackingLeapFrogStraight
+ INTERFACE
+ $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
+ $<INSTALL_INTERFACE:include/include>
+ )
+
+install (FILES ${MODEL_HEADERS} DESTINATION include/${MODEL_NAMESPACE})
+
+# #-- -- -- -- -- -- -- -- -- --
+# #code unit testing
+CORSIKA_ADD_TEST (testTrackingLeapFrogStraight)
+target_link_libraries (
+ testTrackingLeapFrogStraight
+ ProcessTrackingLeapFrogStraight
+ CORSIKAtesting
+)
diff --git a/Processes/TrackingLine/TrackingLine.cc b/Processes/TrackingLeapFrogStraight/Tracking.cc
similarity index 92%
rename from Processes/TrackingLine/TrackingLine.cc
rename to Processes/TrackingLeapFrogStraight/Tracking.cc
index 09c05e2e8..7d2905f8f 100644
--- a/Processes/TrackingLine/TrackingLine.cc
+++ b/Processes/TrackingLeapFrogStraight/Tracking.cc
@@ -13,7 +13,7 @@
#include <corsika/geometry/QuantityVector.h>
#include <corsika/geometry/Sphere.h>
#include <corsika/geometry/Vector.h>
-#include <corsika/process/tracking_line/TrackingLine.h>
+#include <corsika/process/tracking_bfield/Tracking.h>
#include <limits>
#include <stdexcept>
@@ -22,7 +22,7 @@
using namespace corsika::geometry;
using namespace corsika::units::si;
-namespace corsika::process::tracking_line {
+namespace corsika::process::tracking_bfield {
std::optional<std::pair<TimeType, TimeType>> TimeOfIntersection(Line const& line,
Sphere const& sphere) {
@@ -58,4 +58,4 @@ namespace corsika::process::tracking_line {
return n.dot(delta) / c;
}
}
-} // namespace corsika::process::tracking_line
+} // namespace corsika::process::tracking_bfield
diff --git a/Processes/TrackingLeapFrogStraight/Tracking.h b/Processes/TrackingLeapFrogStraight/Tracking.h
new file mode 100644
index 000000000..fd15b89c8
--- /dev/null
+++ b/Processes/TrackingLeapFrogStraight/Tracking.h
@@ -0,0 +1,183 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#pragma once
+
+#include <corsika/geometry/Line.h>
+#include <corsika/geometry/Plane.h>
+#include <corsika/geometry/Sphere.h>
+#include <corsika/geometry/Trajectory.h>
+#include <corsika/geometry/Vector.h>
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/units/PhysicalUnits.h>
+#include <corsika/process/tracking_line/Tracking.h>
+#include <corsika/utl/quartic.h>
+
+#include <type_traits>
+#include <utility>
+
+#include <fstream>
+
+namespace corsika::process {
+
+ namespace tracking_leapfrog_straight {
+
+ typedef corsika::geometry::Vector<corsika::units::si::magnetic_flux_density_d>
+ MagneticFieldVector;
+
+ /**
+ * \class Tracking
+ *
+ * The class tracking_leapfrog_straight::Tracking inherits from
+ * tracking_line::Tracking and adds a (two-step) Leap-Frog
+ * algorithms with two halve-steps and magnetic deflection.
+ *
+ * The two halve steps are implemented as two
+ * `tracking_line::Tracking`s and all geometry intersections are,
+ * thus, based on those two straight line elements.
+ *
+ * As a precaution for numerical instability, the steplength is
+ * limited to correspond to a straight line distance to the next
+ * volume intersection. In typical situations this leads to about
+ * one full leap-frog step to the next volume boundary.
+ *
+ **/
+
+ class Tracking : public tracking_line::Tracking {
+
+ public:
+ template <typename Particle>
+ auto GetTrack(Particle& particle) {
+ using namespace corsika::units::si;
+ using namespace corsika::geometry;
+ geometry::Vector<SpeedType::dimension_type> initialVelocity =
+ particle.GetMomentum() / particle.GetEnergy() * corsika::units::constants::c;
+
+ const Point initialPosition = particle.GetPosition();
+ C8LOG_DEBUG(
+ "TrackingB pid: {}"
+ " , E = {} GeV",
+ particle.GetPID(), particle.GetEnergy() / 1_GeV);
+ C8LOG_DEBUG("TrackingB pos: {}", initialPosition.GetCoordinates());
+ C8LOG_DEBUG("TrackingB E: {} GeV", particle.GetEnergy() / 1_GeV);
+ C8LOG_DEBUG("TrackingB p: {} GeV",
+ particle.GetMomentum().GetComponents() / 1_GeV);
+ C8LOG_DEBUG("TrackingB v: {} ", initialVelocity.GetComponents());
+
+ typedef decltype(particle.GetNode()) node_type;
+ const node_type volumeNode = particle.GetNode();
+ auto magneticfield =
+ volumeNode->GetModelProperties().GetMagneticField(initialPosition);
+
+ // charge of the particle
+ const int chargeNumber = particle.GetChargeNumber();
+ const auto magnitudeB = magneticfield.GetNorm();
+ C8LOG_DEBUG("field={} uT, chargeNumber={}, magnitudeB={} uT",
+ magneticfield.GetComponents() / 1_uT, chargeNumber, magnitudeB / 1_T);
+
+ // we need to limit maximum step-length since we absolutely
+ // need to follow strongly curved trajectories segment-wise,
+ // at least if we don't employ concepts as "Helix
+ // Trajectories" or similar
+ auto const momentumVerticalMag =
+ particle.GetMomentum() -
+ particle.GetMomentum().parallelProjectionOnto(magneticfield);
+ LengthType const gyroradius =
+ (chargeNumber == 0 || magnitudeB == 0_T
+ ? std::numeric_limits<TimeType::value_type>::infinity() * 1_m
+ : momentumVerticalMag.norm() * 1_V /
+ (corsika::units::constants::c * abs(chargeNumber) * magnitudeB *
+ 1_eV));
+ const double maxRadians = 0.01;
+ const LengthType steplimit = 2 * cos(maxRadians) * sin(maxRadians) * gyroradius;
+ C8LOG_DEBUG("gyroradius {}, Steplimit: {}", gyroradius, steplimit);
+
+ // calculate first halve step for "steplimit"
+ const auto initialMomentum = particle.GetMomentum();
+ const auto absMomentum = initialMomentum.norm();
+ const auto absVelocity = initialVelocity.norm();
+ const geometry::Vector<dimensionless_d> direction = initialVelocity.normalized();
+ ;
+ // check if particle is moving at all
+ if (absVelocity * 1_s == 0_m) {
+ return std::make_tuple(
+ geometry::LineTrajectory(geometry::Line(initialPosition, initialVelocity),
+ 0_s),
+ volumeNode);
+ }
+
+ // check, where the first halve-step direction has geometric intersections
+ const auto [initialTrack, initialTrackNextVolume] =
+ tracking_line::Tracking::GetTrack(particle);
+ { [[maybe_unused]] auto& initialTrackNextVolume_dum = initialTrackNextVolume; }
+ const auto initialTrackLength = initialTrack.GetLength(1);
+
+ C8LOG_DEBUG("initialTrack(0)={}, initialTrack(1)={}, initialTrackLength={}",
+ initialTrack.GetPosition(0).GetCoordinates(),
+ initialTrack.GetPosition(1).GetCoordinates(), initialTrackLength);
+
+ // avoid any intersections within first halve steplength
+ LengthType firstHalveSteplength = std::min(steplimit, initialTrackLength) / 2;
+
+ C8LOG_DEBUG("first halve step length {}, steplimit={}, initialTrackLength={}",
+ firstHalveSteplength, steplimit, initialTrackLength);
+ // perform the first halve-step
+ const Point position_mid = initialPosition + direction * firstHalveSteplength;
+ const auto k = chargeNumber * corsika::units::constants::c * 1_eV /
+ (particle.GetMomentum().norm() * 1_V);
+ const auto new_direction =
+ direction + direction.cross(magneticfield) * firstHalveSteplength * 2 * k;
+ const auto new_direction_norm = new_direction.norm(); // by design this is >1
+ C8LOG_DEBUG("position_mid={}, new_direction={}, new_direction_norm={}",
+ position_mid.GetCoordinates(), new_direction.GetComponents(),
+ new_direction_norm);
+
+ // check, where the second halve-step direction has geometric intersections
+ particle.SetPosition(position_mid);
+ particle.SetMomentum(new_direction * absMomentum);
+ const auto [finalTrack, finalTrackNextVolume] =
+ tracking_line::Tracking::GetTrack(particle);
+ particle.SetPosition(initialPosition); // this is not nice...
+ particle.SetMomentum(initialMomentum); // this is not nice...
+ const auto finalTrackLength = finalTrack.GetLength(1);
+
+ C8LOG_DEBUG("finalTrack(0)={}, finalTrack(1)={}, finalTrackLength={}",
+ finalTrack.GetPosition(0).GetCoordinates(),
+ finalTrack.GetPosition(1).GetCoordinates(), finalTrackLength);
+
+ const LengthType secondLeapFrogLength = firstHalveSteplength * new_direction_norm;
+ const LengthType secondHalveStepLength =
+ std::min(secondLeapFrogLength, finalTrackLength);
+
+ // perform the second halve-step
+ const Point finalPosition = position_mid + new_direction * secondHalveStepLength;
+
+ const LengthType totalStep = firstHalveSteplength + secondHalveStepLength;
+ const auto delta_pos = finalPosition - initialPosition;
+ const auto distance = delta_pos.norm();
+
+ return std::make_tuple(
+ geometry::LineTrajectory(
+ geometry::Line(initialPosition,
+ (distance == 0_m ? initialVelocity
+ : delta_pos.normalized() * absVelocity)),
+ distance / absVelocity, // straight distance
+ totalStep / absVelocity, // bend distance
+ initialVelocity,
+ new_direction.normalized() * absVelocity), // trajectory
+ (finalTrackLength > secondLeapFrogLength
+ ? volumeNode
+ : finalTrackNextVolume)); // next step volume
+ }
+ };
+
+ } // namespace tracking_leapfrog_straight
+
+} // namespace corsika::process
diff --git a/Processes/TrackingLeapFrogStraight/testTrackingBFieldStack.h b/Processes/TrackingLeapFrogStraight/testTrackingBFieldStack.h
new file mode 100644
index 000000000..98d7ce3fe
--- /dev/null
+++ b/Processes/TrackingLeapFrogStraight/testTrackingBFieldStack.h
@@ -0,0 +1,58 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#pragma once
+
+#include <corsika/environment/Environment.h>
+
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/Vector.h>
+
+#include <corsika/particles/ParticleProperties.h>
+
+#include <corsika/stack/CombinedStack.h>
+#include <corsika/stack/node/GeometryNodeStackExtension.h>
+#include <corsika/stack/nuclear_extension/NuclearStackExtension.h>
+
+#include <corsika/units/PhysicalUnits.h>
+
+class TestMagneticField {
+ using MagneticFieldVector =
+ corsika::geometry::Vector<corsika::units::si::magnetic_flux_density_d>;
+
+ TestMagneticField() = delete;
+
+public:
+ TestMagneticField(const corsika::units::si::MagneticFluxType& field)
+ : Bz_(field) {}
+
+ void SetMagneticField(const corsika::units::si::MagneticFluxType& field) { Bz_ = field; }
+ MagneticFieldVector GetMagneticField(corsika::geometry::Point const& p) const {
+ using namespace corsika::units::si;
+ return MagneticFieldVector(p.GetCoordinateSystem(), 0_T, 0_T, Bz_);
+ }
+
+private:
+ corsika::units::si::MagneticFluxType Bz_;
+};
+
+using TestEnvironmentType = corsika::environment::Environment<TestMagneticField>;
+
+template <typename T>
+using SetupGeometryDataInterface =
+ corsika::stack::node::GeometryDataInterface<T, TestEnvironmentType>;
+
+// combine particle data stack with geometry information for tracking
+template <typename StackIter>
+using StackWithGeometryInterface = corsika::stack::CombinedParticleInterface<
+ corsika::stack::nuclear_extension::ParticleDataStack::MPIType,
+ SetupGeometryDataInterface, StackIter>;
+
+using TestTrackingBFieldStack = corsika::stack::CombinedStack<
+ typename corsika::stack::nuclear_extension::ParticleDataStack::StackImpl,
+ corsika::stack::node::GeometryData<TestEnvironmentType>, StackWithGeometryInterface>;
diff --git a/Processes/TrackingLeapFrogStraight/testTrackingLeapFrogStraight.cc b/Processes/TrackingLeapFrogStraight/testTrackingLeapFrogStraight.cc
new file mode 100644
index 000000000..e75cbf741
--- /dev/null
+++ b/Processes/TrackingLeapFrogStraight/testTrackingLeapFrogStraight.cc
@@ -0,0 +1,127 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#include <corsika/process/tracking_bfield/Tracking.h>
+
+#include <corsika/particles/ParticleProperties.h>
+
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/Sphere.h>
+#include <corsika/geometry/Vector.h>
+
+#include <corsika/setup/SetupEnvironment.h>
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+using corsika::setup::Trajectory;
+
+#include <catch2/catch.hpp>
+
+using namespace corsika;
+using namespace corsika::particles;
+using namespace corsika::process;
+using namespace corsika::units;
+using namespace corsika::geometry;
+using namespace corsika::units::si;
+
+template <typename T>
+int sgn(T val) {
+ return (T(0) < val) - (val < T(0));
+}
+
+TEST_CASE("TrackingBField") {
+
+ logging::SetLevel(logging::level::trace);
+
+ const HEPEnergyType P0 = 10_GeV;
+
+ auto PID = GENERATE(as<Code>{}, Code::MuPlus, Code::MuPlus, Code::Gamma);
+ auto Bfield = GENERATE(as<MagneticFluxType>{}, 0_T, 50_uT, -50_uT);
+ auto outer = GENERATE(as<bool>{}, true, false);
+
+ SECTION(fmt::format("Tracking PID={}, Bfield={} uT, from outside={}", PID,
+ Bfield / 1_uT, outer)) {
+
+ C8LOG_DEBUG(
+ "********************\n TEST section PID={}, Bfield={} "
+ "uT, outer (?)={}",
+ PID, Bfield / 1_uT, outer);
+
+ const int chargeNumber = GetChargeNumber(PID);
+ LengthType radius = 10_m;
+ int deflect = 0;
+ if (chargeNumber != 0 and Bfield != 0_T) {
+ deflect = -sgn(chargeNumber) * sgn(Bfield / 1_T); // direction of deflection
+ LengthType const gyroradius =
+ P0 * 1_V / (constants::c * abs(chargeNumber) * abs(Bfield) * 1_eV);
+ radius = gyroradius;
+ }
+
+ auto [env, csPtr, worldPtr] = setup::testing::setupEnvironment(Code::Oxygen, Bfield);
+ { [[maybe_unused]] const auto& env_dummy = env; }
+ auto const& cs = *csPtr;
+
+ tracking_leapfrog_straight::Tracking tracking;
+ using tracking_leapfrog_straight::MagneticFieldVector;
+ Point const center(cs, {0_m, 0_m, 0_m});
+ auto target = setup::Environment::CreateNode<geometry::Sphere>(center, radius);
+
+ using MyHomogeneousModel =
+ environment::MediumPropertyModel<environment::UniformMagneticField<
+ environment::HomogeneousMedium<setup::EnvironmentInterface>>>;
+
+ MagneticFieldVector magneticfield(cs, 0_T, 0_T, Bfield);
+ target->SetModelProperties<MyHomogeneousModel>(
+ environment::Medium::AirDry1Atm, magneticfield, 1_g / (1_m * 1_m * 1_m),
+ environment::NuclearComposition(std::vector<Code>{Code::Oxygen},
+ std::vector<float>{1.}));
+ auto* targetPtr = target.get();
+ worldPtr->AddChild(std::move(target));
+
+ auto [stack, viewPtr] = setup::testing::setupStack(PID, 0, 0, P0, targetPtr, cs);
+ { [[maybe_unused]] auto& viewPtr_dum = viewPtr; }
+ auto particle = stack->first();
+ // Note: momentum in X-direction
+ // magnetic field in Z-direction
+ // put particle on x_start, 0, 0
+ // expect intersections somewere in +-y_start
+
+ if (outer) {
+ particle.SetNode(worldPtr); // set particle inside "target" volume
+ } else {
+ particle.SetNode(targetPtr); // set particle outside "target" volume
+ }
+ particle.SetPosition(Point(cs, -radius, 0_m, 0_m));
+
+ auto [traj, nextVol] = tracking.GetTrack(particle);
+ particle.SetNode(nextVol);
+ particle.SetPosition(traj.GetPosition(1));
+ particle.SetMomentum(traj.GetDirection(1) * particle.GetMomentum().norm());
+ if (outer) {
+ // now we know we are in target volume, depending on "outer"
+ CHECK(traj.GetLength(1) == 0_m);
+ CHECK(nextVol == targetPtr);
+ }
+ // move forward, until we leave target volume
+ while (nextVol == targetPtr) {
+ const auto [traj2, nextVol2] = tracking.GetTrack(particle);
+ nextVol = nextVol2;
+ particle.SetNode(nextVol);
+ particle.SetPosition(traj2.GetPosition(1));
+ particle.SetMomentum(traj2.GetDirection(1) * particle.GetMomentum().norm());
+ }
+ CHECK(nextVol == worldPtr);
+
+ Point pointCheck(cs, (deflect == 0 ? radius : 0_m), (deflect * radius), 0_m);
+
+ C8LOG_DEBUG("testTrackingLineStack: deflect={}, momentum={}, pos={}, pos_check={}",
+ deflect, particle.GetMomentum().GetComponents(),
+ particle.GetPosition().GetCoordinates(), pointCheck.GetCoordinates());
+
+ CHECK((particle.GetPosition() - pointCheck).norm() / radius == Approx(0).margin(1e-3));
+ }
+}
diff --git a/Processes/TrackingLine/CMakeLists.txt b/Processes/TrackingLine/CMakeLists.txt
index 36d53e30c..a5376eac1 100644
--- a/Processes/TrackingLine/CMakeLists.txt
+++ b/Processes/TrackingLine/CMakeLists.txt
@@ -1,12 +1,6 @@
set (
MODEL_HEADERS
- TrackingLine.h
- dump_bh.hpp
- )
-
-set (
- MODEL_SOURCES
- TrackingLine.cc
+ Tracking.h
)
set (
@@ -14,19 +8,14 @@ set (
corsika/process/tracking_line
)
-add_library (ProcessTrackingLine STATIC ${MODEL_SOURCES})
+add_library (ProcessTrackingLine INTERFACE)
CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessTrackingLine ${MODEL_NAMESPACE} ${MODEL_HEADERS})
-set_target_properties (
- ProcessTrackingLine
- PROPERTIES
- VERSION ${PROJECT_VERSION}
- SOVERSION 1
- )
-
# target dependencies on other libraries (also the header onlys)
target_link_libraries (
ProcessTrackingLine
+ INTERFACE
+ ProcessTrackingIntersects
CORSIKAsetup
CORSIKAutilities
CORSIKAenvironment
diff --git a/Processes/TrackingLine/Tracking.cc b/Processes/TrackingLine/Tracking.cc
new file mode 100644
index 000000000..01cd6c134
--- /dev/null
+++ b/Processes/TrackingLine/Tracking.cc
@@ -0,0 +1,26 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#include <corsika/environment/Environment.h>
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/QuantityVector.h>
+#include <corsika/geometry/Sphere.h>
+#include <corsika/geometry/Vector.h>
+#include <corsika/process/tracking_line/Tracking.h>
+#include <corsika/logging/Logging.h>
+
+#include <limits>
+#include <stdexcept>
+#include <utility>
+
+using namespace corsika::geometry;
+using namespace corsika::units::si;
+
+namespace corsika::process::tracking_line {
+
+} // namespace corsika::process::tracking_line
diff --git a/Processes/TrackingLine/Tracking.h b/Processes/TrackingLine/Tracking.h
new file mode 100644
index 000000000..799f2194f
--- /dev/null
+++ b/Processes/TrackingLine/Tracking.h
@@ -0,0 +1,126 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#pragma once
+
+#include <corsika/geometry/Line.h>
+#include <corsika/geometry/Plane.h>
+#include <corsika/geometry/Sphere.h>
+#include <corsika/geometry/Trajectory.h>
+#include <corsika/geometry/Vector.h>
+#include <corsika/geometry/Intersections.hpp>
+#include <corsika/units/PhysicalUnits.h>
+#include <corsika/logging/Logging.h>
+#include <corsika/setup/SetupEnvironment.h>
+#include <corsika/process/tracking/Intersect.hpp>
+
+#include <type_traits>
+#include <utility>
+
+namespace corsika::process {
+
+ namespace tracking_line {
+
+ /**
+ * \class Tracking
+ *
+ *
+ *
+ **/
+
+ class Tracking : public corsika::process::tracking::Intersect<Tracking> {
+
+ public:
+ template <typename TParticle>
+ auto GetTrack(TParticle const& particle) {
+ using namespace corsika::units::si;
+ using namespace corsika::geometry;
+ geometry::Vector<SpeedType::dimension_type> const initialVelocity =
+ particle.GetMomentum() / particle.GetEnergy() * corsika::units::constants::c;
+
+ auto const initialPosition = particle.GetPosition();
+ C8LOG_DEBUG(
+ "Tracking pid: {}"
+ " , E = {} GeV",
+ particle.GetPID(), particle.GetEnergy() / 1_GeV);
+ C8LOG_DEBUG("Tracking pos: {}", initialPosition.GetCoordinates());
+ C8LOG_DEBUG("Tracking E: {} GeV", particle.GetEnergy() / 1_GeV);
+ C8LOG_DEBUG("Tracking p: {} GeV",
+ particle.GetMomentum().GetComponents() / 1_GeV);
+ C8LOG_DEBUG("Tracking v: {} ", initialVelocity.GetComponents());
+
+ // traverse the environment volume tree and find next
+ // intersection
+ auto [minTime, minNode] = tracking::Intersect<Tracking>::nextIntersect(particle);
+
+ return std::make_tuple(
+ geometry::LineTrajectory(geometry::Line(initialPosition, initialVelocity),
+ minTime), // trajectory
+ minNode); // next volume node
+ }
+
+
+
+ template <typename TParticle, typename TMedium>
+ static geometry::Intersections Intersect(const TParticle& particle,
+ const corsika::geometry::Sphere& sphere,
+ const TMedium&) {
+ using namespace corsika::units::si;
+ auto const delta = particle.GetPosition() - sphere.GetCenter();
+ auto const velocity =
+ particle.GetMomentum() / particle.GetEnergy() * corsika::units::constants::c;
+ auto const vSqNorm = velocity.squaredNorm();
+ auto const R = sphere.GetRadius();
+
+ auto const vDotDelta = velocity.dot(delta);
+ auto const discriminant =
+ vDotDelta * vDotDelta - vSqNorm * (delta.squaredNorm() - R * R);
+
+ if (discriminant.magnitude() > 0) {
+ auto const sqDisc = sqrt(discriminant);
+ auto const invDenom = 1 / vSqNorm;
+ return geometry::Intersections((-vDotDelta - sqDisc) * invDenom,
+ (-vDotDelta + sqDisc) * invDenom);
+ }
+ return geometry::Intersections();
+ }
+
+ template <typename TParticle, typename TBaseNodeType>
+ static geometry::Intersections Intersect(const TParticle& particle,
+ const TBaseNodeType& volumeNode) {
+ const geometry::Sphere* sphere =
+ dynamic_cast<const geometry::Sphere*>(&volumeNode.GetVolume());
+ if (sphere) {
+ return Intersect(particle, *sphere, volumeNode.GetModelProperties());
+ }
+ throw std::runtime_error(
+ "The Volume type provided is not supported in Intersect(particle, node)");
+ }
+
+
+ template <typename TParticle, typename TMedium>
+ static geometry::Intersections Intersect(const TParticle& particle,
+ const geometry::Plane& plane,
+ const TMedium& medium) {
+ using namespace corsika::units::si;
+ auto const delta = plane.GetCenter() - particle.GetPosition();
+ auto const velocity =
+ particle.GetMomentum() / particle.GetEnergy() * corsika::units::constants::c;
+ auto const n = plane.GetNormal();
+ auto const c = n.dot(velocity);
+
+ return Intersections(
+ c.magnitude() == 0 ? std::numeric_limits<TimeType::value_type>::infinity() * 1_s
+ : n.dot(delta) / c);
+ }
+
+ };
+
+ } // namespace tracking_line
+
+} // namespace corsika::process
diff --git a/Processes/TrackingLine/TrackingLine.h b/Processes/TrackingLine/TrackingLine.h
deleted file mode 100644
index e7e98cb53..000000000
--- a/Processes/TrackingLine/TrackingLine.h
+++ /dev/null
@@ -1,694 +0,0 @@
-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#ifndef _include_corsika_processes_TrackingLine_h_
-#define _include_corsika_processes_TrackingLine_h_
-
-#include <corsika/geometry/Line.h>
-#include <corsika/geometry/Plane.h>
-#include <corsika/geometry/Sphere.h>
-#include <corsika/geometry/Trajectory.h>
-#include <corsika/geometry/Vector.h>
-#include <corsika/particles/ParticleProperties.h>
-#include <corsika/units/PhysicalUnits.h>
-#include <corsika/utl/quartic.h>
-#include <optional>
-#include <type_traits>
-#include <utility>
-
-#include <fstream>
-#include <iostream>
-#include <boost/histogram.hpp>
-#include <boost/histogram/ostream.hpp>
-#include <corsika/process/tracking_line/dump_bh.hpp>
-
-namespace corsika::environment {
- template <typename IEnvironmentModel>
- class Environment;
- template <typename IEnvironmentModel>
- class VolumeTreeNode;
-} // namespace corsika::environment
-
-using namespace boost::histogram;
-/*static auto histL = make_histogram(axis::regular<>(100, 0, 100000, "Leap-Frog-ength L'"));
-static auto histS = make_histogram(axis::regular<>(100, 0, 100000, "Direct Length S"));
-static auto histB = make_histogram(axis::regular<>(100, 0, 100000, "Arc Length B")); */
-static auto histLlog = make_histogram(axis::regular<double, axis::transform::log>(100, 1e-3, 1e7, "Leap-Frog-ength L'"));
-static auto histLloggeo = make_histogram(axis::regular<double, axis::transform::log>(100, 1e-3, 1e7, "Leap-Frog-length L for geometric steps"));
-static auto histLlogmag = make_histogram(axis::regular<double, axis::transform::log>(100, 1e-3, 1e7, "Leap-Frog-length L for magnetic steps"));
-static auto histSlog = make_histogram(axis::regular<double, axis::transform::log>(100, 1e-3, 1e7, "Direct Length S"));
-static auto histBlog = make_histogram(axis::regular<double, axis::transform::log>(100, 1e-3, 1e7, "Arc Length B"));
-/*static auto histLB = make_histogram(axis::regular<>(100, 0, 100, "L - B"));
-static auto histLS = make_histogram(axis::regular<>(100, 0, 100, "L - S"));*/
-static auto histLBrelgeo = make_histogram(axis::regular<double, axis::transform::log> (40,1e-12,1e-2,"L/B -1"));
-static auto histLBrelmag = make_histogram(axis::regular<double, axis::transform::log> (40,1e-12,1e-2,"L/B -1"));
-static auto histLSrelgeo = make_histogram(axis::regular<double, axis::transform::log>(40,1e-12,1e-2, "L/S -1"));
-static auto histLSrelmag = make_histogram(axis::regular<double, axis::transform::log>(40,1e-12,1e-2, "L/S -1"));
-static auto histELSrel = make_histogram(axis::regular<double, axis::transform::log>(50,1e-8,1e-3, "L/S -1"),axis::regular<double, axis::transform::log>(30, 3, 3e1, "E / GeV"));
-static auto histELSrelp = make_histogram(axis::regular<double, axis::transform::log>(50,1e-8,1e-3, "L/S -1"),axis::regular<double, axis::transform::log>(30, 3, 3e1, "E / GeV"));
-static auto histELSrelpi = make_histogram(axis::regular<double, axis::transform::log>(50,1e-8,1e-3, "L/S -1"),axis::regular<double, axis::transform::log>(30, 3, 3e1, "E / GeV"));
-static auto histELSrelmu = make_histogram(axis::regular<double, axis::transform::log>(50,1e-8,1e-3, "L/S -1"),axis::regular<double, axis::transform::log>(30, 3, 3e1, "E / GeV"));
-static auto histELSrele = make_histogram(axis::regular<double, axis::transform::log>(50,1e-8,1e-3, "L/S -1"),axis::regular<double, axis::transform::log>(30, 3, 3e1, "E / GeV"));
-//static auto histBS = make_histogram(axis::regular<>(100, 0, 0.01, "B - S"));
-static auto histLpgeo = make_histogram(axis::regular<double, axis::transform::log>(100, 1e-3, 1e7, "L' für Protonen"));
-static auto histLpigeo = make_histogram(axis::regular<double, axis::transform::log>(100, 1e-3, 1e7, "L' für Pionen"));
-static auto histLmugeo = make_histogram(axis::regular<double, axis::transform::log>(100, 1e-3, 1e7, "L' für Myonen"));
-static auto histLegeo = make_histogram(axis::regular<double, axis::transform::log>(100, 1e-3, 1e7, "L' für Elektronen"));
-static auto histLygeo = make_histogram(axis::regular<double, axis::transform::log>(100, 1e-3, 1e7, "L' für Photonen"));
-static auto histLpmag = make_histogram(axis::regular<double, axis::transform::log>(100, 1e-3, 1e7, "L' für Protonen"));
-static auto histLpimag = make_histogram(axis::regular<double, axis::transform::log>(100, 1e-3, 1e7, "L' für Pionen"));
-static auto histLmumag = make_histogram(axis::regular<double, axis::transform::log>(100, 1e-3, 1e7, "L' für Myonen"));
-static auto histLemag = make_histogram(axis::regular<double, axis::transform::log>(100, 1e-3, 1e7, "L' für Elektronen"));
-static auto histLymag = make_histogram(axis::regular<double, axis::transform::log>(100, 1e-3, 1e7, "L' für Photonen"));
-static double L = 0;
-static std::vector<double> Lsmall;
-static std::vector<double> Bbig;
-static std::vector<double> Sbig;
-static std::vector<double> Lsmallmag;
-static std::vector<double> Bbigmag;
-static std::vector<double> Sbigmag;
-static std::vector<double> Emag;
-static std::vector<double> E;
-static std::vector<double> Steplimitmag;
-
-
-namespace corsika::process {
-
- namespace tracking_line {
-
- std::optional<std::pair<corsika::units::si::TimeType, corsika::units::si::TimeType>>
- TimeOfIntersection(geometry::Line const&, geometry::Sphere const&);
-
- corsika::units::si::TimeType TimeOfIntersection(geometry::Line const&,
- geometry::Plane const&);
-
- class TrackingLine {
-
- public:
- TrackingLine() {};
- ~TrackingLine() {
- /*std::ofstream myfile;
- myfile.open ("histograms.txt");
- myfile << histLlog << std::endl;
- myfile << histBlog << std::endl;
- myfile << histLB << std::endl;
- myfile << histLS << std::endl;
- myfile.close();
-
- std::ofstream myfile2;
- myfile2.open ("lengen.txt");
- myfile2 << "L " << std::endl;
- for(double n : Lsmall) {
- myfile2 << n << '\n';
- }
- myfile2 << "B " << std::endl;
- for(double n : Bbig) {
- myfile2 << n << '\n';
- }
- myfile2 << "S " << std::endl;
- for(double n : Sbig) {
- myfile2 << n << '\n';
- }
- myfile2 << "E in GeV" << std::endl;
- for(double n : E) {
- myfile2 << n << '\n';
- }
- myfile2 << "L mag" << std::endl;
- for(double n : Lsmallmag) {
- myfile2 << n << '\n';
- }
- myfile2 << "B mag" << std::endl;
- for(double n : Bbigmag) {
- myfile2 << n << '\n';
- }
- myfile2 << "S mag" << std::endl;
- for(double n : Sbigmag) {
- myfile2 << n << '\n';
- }
- myfile2 << "E mag in GeV" << std::endl;
- for(double n : Emag) {
- myfile2 << n << '\n';
- }
- myfile2 << "Schrittlänge" << std::endl;
- for(double n : Steplimitmag) {
- myfile2 << n << '\n';
- }
- myfile2.close();
-
- std::ofstream file1("histL.json");
- dump_bh(file1, histL);
- file1.close();
- std::ofstream file2("histS.json");
- dump_bh(file2, histS);
- file2.close();
- std::ofstream file3("histB.json");
- dump_bh(file3, histB);
- file3.close();
- std::ofstream file4("histLB.json");
- dump_bh(file4, histLB);
- file4.close();
- std::ofstream file5("histLS.json");
- dump_bh(file5, histLS);
- file5.close();
- std::ofstream file6("histBS.json");
- dump_bh(file6, histBS);
- file6.close();*/
- /*std::ofstream file7("histLBrelgeo.json");
- dump_bh(file7, histLBrelgeo);
- file7.close();
- std::ofstream file8("histLSrelgeo.json");
- dump_bh(file8, histLSrelgeo);
- file8.close();
- std::ofstream file9("histLBrelmag.json");
- dump_bh(file9, histLBrelmag);
- file9.close();
- std::ofstream file0("histLSrelmag.json");
- dump_bh(file0, histLSrelmag);
- file0.close();
-
- std::ofstream file10("histELSrel.json");
- dump_bh(file10, histELSrel);
- file10.close();
- std::ofstream file11("histLmugeo.json");
- dump_bh(file11, histLmugeo);
- file11.close();
- std::ofstream file12("histLpigeo.json");
- dump_bh(file12, histLpigeo);
- file12.close();
- std::ofstream file13("histLpgeo.json");
- dump_bh(file13, histLpgeo);
- file13.close();
- std::ofstream file14("histLlog.json");
- dump_bh(file14, histLlog);
- file14.close();
- std::ofstream file15("histBlog.json");
- dump_bh(file15, histBlog);
- file15.close();
- std::ofstream file16("histSlog.json");
- dump_bh(file16, histSlog);
- file16.close();
- std::ofstream file17("histELSrelmu.json");
- dump_bh(file17, histELSrelmu);
- file17.close();
- std::ofstream file18("histELSrelp.json");
- dump_bh(file18, histELSrelp);
- file18.close();
- std::ofstream file19("histELSrelpi.json");
- dump_bh(file19, histELSrelpi);
- file19.close();
-
- std::ofstream file21("histLgeo.json");
- dump_bh(file21, histLloggeo);
- file21.close();
- std::ofstream file22("histLmag.json");
- dump_bh(file22, histLlogmag);
- file22.close();
- std::ofstream file23("histLmumag.json");
- dump_bh(file23, histLmumag);
- file23.close();
- std::ofstream file24("histLpimag.json");
- dump_bh(file24, histLpimag);
- file24.close();
- std::ofstream file25("histLpmag.json");
- dump_bh(file25, histLpmag);
- file25.close();
- std::ofstream file26("histLemag.json");
- dump_bh(file26, histLemag);
- file26.close();
- std::ofstream file27("histLegeo.json");
- dump_bh(file27, histLegeo);
- file27.close();
- std::ofstream file28("histELSrele.json");
- dump_bh(file28, histELSrele);
- file28.close();
- std::ofstream file29("histLymag.json");
- dump_bh(file29, histLymag);
- file29.close();
- std::ofstream file20("histLygeo.json");
- dump_bh(file20, histLygeo);
- file20.close();*/
-
- };
-
- template <typename Particle> // was Stack previously, and argument was
- // Stack::StackIterator
- auto GetTrack(Particle const& p) {
- using namespace corsika::units::si;
- using namespace corsika::geometry;
- geometry::Vector<SpeedType::dimension_type> velocity =
- p.GetMomentum() / p.GetEnergy() * corsika::units::constants::c;
-
- auto const currentPosition = p.GetPosition();
- std::cout << "TrackingLine pid: " << p.GetPID()
- << " , E = " << p.GetEnergy() / 1_GeV << " GeV" << std::endl;
- std::cout << "TrackingLine pos: "
- << currentPosition.GetCoordinates()
- // << " [" << p.GetNode()->GetModelProperties().GetName() << "]"
- << std::endl;
- std::cout << "TrackingLine E: " << p.GetEnergy() / 1_GeV << " GeV" << std::endl;
- std::cout << "TrackingLine p: " << p.GetMomentum().GetComponents() / 1_GeV
- << " GeV " << std::endl;
- std::cout << "TrackingLine v: " << velocity.GetComponents() << std::endl;
-
- geometry::Line line(currentPosition, velocity);
-
- auto const* currentLogicalVolumeNode = p.GetNode();
- //~ auto const* currentNumericalVolumeNode =
- //~ fEnvironment.GetUniverse()->GetContainingNode(currentPosition);
- auto const numericallyInside =
- currentLogicalVolumeNode->GetVolume().Contains(currentPosition);
-
- std::cout << "numericallyInside = " << (numericallyInside ? "true" : "false") << std::endl;
-
- auto const& children = currentLogicalVolumeNode->GetChildNodes();
- auto const& excluded = currentLogicalVolumeNode->GetExcludedNodes();
-
- std::vector<std::pair<TimeType, decltype(p.GetNode())>> intersections;
-
- //charge of the particle
- int chargeNumber = 0;
- if (corsika::particles::IsNucleus(p.GetPID())) {
- chargeNumber = p.GetNuclearZ();
- }
- else {
- chargeNumber = corsika::particles::GetChargeNumber(p.GetPID());
- }
- auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(currentPosition);
- std::cout << " Magnetic Field: " << magneticfield.GetComponents() / 1_uT << " uT " << std::endl;
- LengthType Steplength1 = 0_m;
- LengthType Steplength2 = 0_m;
- LengthType Steplimit = 1_m / 0;
- //infinite kann man auch anders schreiben
- bool intersection = true;
-
- // for entering from outside
- auto addIfIntersects = [&](auto const& vtn) {
- auto const& volume = vtn.GetVolume();
- auto const& sphere = dynamic_cast<geometry::Sphere const&>(
- volume); // for the moment we are a bit bold here and assume
- // everything is a sphere, crashes with exception if not
-
- // creating Line with magnetic field
- if (chargeNumber != 0) {
- auto k = chargeNumber * corsika::units::constants::c * 1_eV / (p.GetMomentum().norm() * 1_V);
- geometry::Vector<dimensionless_d> const directionBefore = velocity.normalized();
- // determine steplength to next volume
- double a = ((directionBefore.cross(magneticfield)).dot(currentPosition - sphere.GetCenter()) * k + 1) * 4 /
- (1_m * 1_m * (directionBefore.cross(magneticfield)).GetSquaredNorm() * k * k);
- double b = directionBefore.dot(currentPosition - sphere.GetCenter()) * 8 /
- ((directionBefore.cross(magneticfield)).GetSquaredNorm() * k * k * 1_m * 1_m * 1_m);
- double c = ((currentPosition - sphere.GetCenter()).GetSquaredNorm() -
- (sphere.GetRadius() * sphere.GetRadius())) * 4 /
- ((directionBefore.cross(magneticfield)).GetSquaredNorm() * k * k * 1_m * 1_m * 1_m * 1_m);
- std::complex<double>* solutions = solve_quartic(0, a, b, c);
- std::vector<double> tmp;
- for (int i = 0; i < 4; i++) {
- if (solutions[i].imag() == 0 && solutions[i].real() > 0) {
- tmp.push_back(solutions[i].real());
- std::cout << "Solutions for next Volume: " << solutions[i].real() << std::endl;
- }
- }
- if (tmp.size() > 0) {
- Steplength1 = 1_m * *std::min_element(tmp.begin(), tmp.end());
- std::cout << "Steplength to next volume = " << Steplength1 << std::endl;
- std::cout << "Time to next volume = " << Steplength1 / velocity.norm() << std::endl;
- }
- else {
- std::cout << "no intersection (1)!" << std::endl;
- intersection = false;
- }
- delete[] solutions;
- }
-
- if (intersection) {
- auto line1 = MagneticStep(p, line, Steplength1, false);
- // instead of changing the line with magnetic field, the TimeOfIntersection() could be changed
- // using line has some errors for huge steps
-
- if (auto opt = TimeOfIntersection(line1, sphere); opt.has_value()) {
- auto const[t1, t2] = *opt;
- std::cout << "intersection times: " << t1 / 1_s << "; "
- << t2 / 1_s
- // << " " << vtn.GetModelProperties().GetName()
- << std::endl;
- if (t1.magnitude() > 0)
- intersections.emplace_back(t1, &vtn);
- else if (t2.magnitude() > 0)
- std::cout << "inside other volume" << std::endl;
- }
- }
- };
-
- for (auto const& child : children) { addIfIntersects(*child); }
- for (auto const* ex : excluded) { addIfIntersects(*ex); }
-
- {
- auto const& sphere = dynamic_cast<geometry::Sphere const&>(
- currentLogicalVolumeNode->GetVolume());
- // for the moment we are a bit bold here and assume
- // everything is a sphere, crashes with exception if not
-
- // creating Line with magnetic field
- if (chargeNumber != 0) {
- auto k = chargeNumber * corsika::units::constants::c * 1_eV / (p.GetMomentum().norm() * 1_V);
- geometry::Vector<dimensionless_d> const directionBefore = velocity.normalized();
- // determine steplength to next volume
- double a = ((directionBefore.cross(magneticfield)).dot(currentPosition - sphere.GetCenter()) * k + 1) * 4 /
- (1_m * 1_m * (directionBefore.cross(magneticfield)).GetSquaredNorm() * k * k);
- double b = directionBefore.dot(currentPosition - sphere.GetCenter()) * 8 /
- ((directionBefore.cross(magneticfield)).GetSquaredNorm() * k * k * 1_m * 1_m * 1_m);
- double c = ((currentPosition - sphere.GetCenter()).GetSquaredNorm() -
- (sphere.GetRadius() * sphere.GetRadius())) * 4 /
- ((directionBefore.cross(magneticfield)).GetSquaredNorm() * k * k * 1_m * 1_m * 1_m * 1_m);
- std::complex<double>* solutions = solve_quartic(0, a, b, c);
- std::vector<double> tmp;
- for (int i = 0; i < 4; i++) {
- if (solutions[i].imag() == 0 && solutions[i].real() > 0) {
- tmp.push_back(solutions[i].real());
- std::cout << "Solutions for current Volume: " << solutions[i].real() << std::endl;
- }
- }
- if (tmp.size() > 0) {
- Steplength2 = 1_m * *std::min_element(tmp.begin(), tmp.end());
- if (numericallyInside == false) {
- int p = std::min_element(tmp.begin(), tmp.end()) - tmp.begin();
- tmp.erase(tmp.begin() + p);
- Steplength2 = 1_m * *std::min_element(tmp.begin(), tmp.end());
- }
- std::cout << "steplength out of current volume = " << Steplength2 << std::endl;
- std::cout << "Time out of current volume = " << Steplength2 / velocity.norm() << std::endl;
- }
- else {
- std::cout << "no intersection (2)!" << std::endl;
- // what to do when this happens? (very unlikely)
- }
- delete[] solutions;
-
- auto const momentumVerticalMag = p.GetMomentum() -
- p.GetMomentum().parallelProjectionOnto(magneticfield);
- LengthType const gyroradius = momentumVerticalMag.norm() * 1_V /
- (corsika::units::constants::c * abs(chargeNumber) *
- magneticfield.GetNorm() * 1_eV);
- Steplimit = 2 * sin(0.1) * gyroradius;
- std::cout << "Steplimit: " << Steplimit << std::endl;
- }
-
- auto line2 = MagneticStep(p, line, Steplength2, false);
- // instead of changing the line with magnetic field, the TimeOfIntersection() could be changed
- // using line has some errors for huge steps
-
- [[maybe_unused]] auto const[t1, t2] = *TimeOfIntersection(line2, sphere);
- [[maybe_unused]] auto dummy_t1 = t1;
- intersections.emplace_back(t2, currentLogicalVolumeNode->GetParent());
- }
-
- auto const minIter = std::min_element(
- intersections.cbegin(), intersections.cend(),
- [](auto const& a, auto const& b) { return a.first < b.first; });
-
- TimeType min;
-
- if (minIter == intersections.cend()) {
- min = 1_s; // todo: do sth. more reasonable as soon as tracking is able
- // to handle the numerics properly
- throw std::runtime_error("no intersection with anything!");
- }
- else {
- min = minIter->first;
- }
-
- std::cout << " t-intersect: "
- << min
- // << " " << minIter->second->GetModelProperties().GetName()
- << std::endl;
-
- //if the used Steplengths are too big, min gets false and we do another Step
- std::cout << Steplength1 << " " << Steplength2 << std::endl;
- std::cout << intersection << " " << Steplimit << std::endl;
-
- std::cout << "Test: (1) " << (Steplength1 > Steplimit || Steplength1 < 1e-6_m) << std::endl;
- std::cout << "Test: (2) " << (Steplength2 > Steplimit) << std::endl;
-
- if ((Steplength1 > Steplimit || Steplength1 < 1e-6_m) && Steplength2 > Steplimit) {
- min = Steplimit / velocity.norm();
- auto lineWithB = MagneticStep(p, line, Steplimit, true);
-
- //histL(L);
- histLlog(L);
- histLlogmag(L);
- //histS(lineWithB.ArcLength(0_s,min) / 1_m);
- histSlog(lineWithB.ArcLength(0_s, min) / 1_m);
- //histLS(L-lineWithB.ArcLength(0_s,min) / 1_m);
-
- std::cout << "is it here? TrackingLine.h (2)" << std::endl;
-
- if (chargeNumber != 0) {
- if (L * 1_m / lineWithB.ArcLength(0_s, min) - 1 > 0) {
- histLSrelmag(L * 1_m / lineWithB.ArcLength(0_s, min) - 1);
- histELSrel(L * 1_m / lineWithB.ArcLength(0_s, min) - 1, p.GetEnergy() / 1_GeV);
- }
- else {
- Lsmallmag.push_back(L); Sbigmag.push_back(lineWithB.ArcLength(0_s, min) / 1_m);
- Emag.push_back(p.GetEnergy() / 1_GeV); Steplimitmag.push_back(Steplimit / 1_m);
- }
- auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(currentPosition);
- geometry::Vector<SpeedType::dimension_type> const velocityVerticalMag = velocity -
- velocity.parallelProjectionOnto(magneticfield);
- LengthType const gyroradius = p.GetEnergy() * velocityVerticalMag.norm() * 1_V /
- (corsika::units::constants::cSquared * abs(chargeNumber) *
- magneticfield.GetNorm() * 1_eV);
- std::cout << "Schrittlaenge " << velocity.norm() * min << " gyroradius " << gyroradius << std::endl;
- LengthType B = 2 * gyroradius * asin(lineWithB.ArcLength(0_s, min) / 2 / gyroradius);
- std::cout << "Bogenlaenge" << B << std::endl;
- //histB(B / 1_m);
- histBlog(B / 1_m);
- if (L - B / 1_m > 0) {
- //histLB(L-B/1_m);
- //histBS(B/1_m - lineWithB.ArcLength(0_s,min) / 1_m);
- if(B != 0_m)
- histLBrelmag(L * 1_m / B - 1);
- }
- else {
- Bbigmag.push_back(B / 1_m);
- }
- }
- else {
- //histB(lineWithB.ArcLength(0_s,min) / 1_m);
- histBlog(lineWithB.ArcLength(0_s, min) / 1_m);
- }
-
- int pdg = static_cast<int>(particles::GetPDG(p.GetPID()));
- if (abs(pdg) == 13) {
- histLmumag(L);
- histELSrelmu(L * 1_m / lineWithB.ArcLength(0_s, min) - 1, p.GetEnergy() / 1_GeV);
- }
- if (abs(pdg) == 211 || abs(pdg) == 111) {
- histLpimag(L);
- if (chargeNumber != 0)
- histELSrelpi(L * 1_m / lineWithB.ArcLength(0_s, min) - 1, p.GetEnergy() / 1_GeV);
- }
- if (abs(pdg) == 2212 || abs(pdg) == 2112) {
- histLpmag(L);
- if (chargeNumber != 0)
- histELSrelp(L * 1_m / lineWithB.ArcLength(0_s, min) - 1, p.GetEnergy() / 1_GeV);
- }
- if (abs(pdg) == 11) {
- histLemag(L);
- if (abs(pdg) == 22) {
- histLymag(L);
- }
- }
- return std::make_tuple(geometry::Trajectory<geometry::Line>(line, min),
- Steplimit * 1.0001, Steplimit, minIter->second);
- }
-
- std::cout << "is it here? TrackingLine.h (3)" << std::endl;
-
- auto lineWithB = MagneticStep(p, line, velocity.norm() * min, true);
-
- //histL(L);
- histLlog(L);
- histLloggeo(L);
- //histS(lineWithB.ArcLength(0_s,min) / 1_m);
- histSlog(lineWithB.ArcLength(0_s, min) / 1_m);
- //histLS(L-lineWithB.ArcLength(0_s,min) / 1_m);
-
- if (chargeNumber != 0) {
- if (L * 1_m / lineWithB.ArcLength(0_s, min) - 1 > 0) {
- histLSrelgeo(L * 1_m / lineWithB.ArcLength(0_s, min) - 1);
- histELSrel(L * 1_m / lineWithB.ArcLength(0_s, min) - 1, p.GetEnergy() / 1_GeV);
- }
- else {
- Lsmall.push_back(L); Sbig.push_back(lineWithB.ArcLength(0_s, min) / 1_m);
- E.push_back(p.GetEnergy() / 1_GeV);
- }
- auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(currentPosition);
- geometry::Vector<SpeedType::dimension_type> const velocityVerticalMag = velocity -
- velocity.parallelProjectionOnto(magneticfield);
- LengthType const gyroradius = p.GetEnergy() * velocityVerticalMag.norm() * 1_V /
- (corsika::units::constants::cSquared * abs(chargeNumber) *
- magneticfield.GetNorm() * 1_eV);
- std::cout << "Schrittlaenge " << velocity.norm() * min << " gyroradius " << gyroradius << std::endl;
-
- //irgendetwas geht schief: Schrittlänge viel größer als gyroradius
- if (lineWithB.ArcLength(0_s, min) < gyroradius) {
- LengthType B = 2 * gyroradius * asin(lineWithB.ArcLength(0_s, min) / 2 / gyroradius);
- std::cout << "Bogenlaenge" << B << std::endl;
- //histB(B / 1_m);
- histBlog(B / 1_m);
- if (L - B / 1_m > 0) {
- //histLB(L-B/1_m);
- //histBS(B/1_m - lineWithB.ArcLength(0_s,min) / 1_m);
- if (B != 0_m)
- histLBrelgeo(L * 1_m / B - 1);
- }
- else {
- Bbig.push_back(B / 1_m);
- }
- }
- }
- else {
- //histB(lineWithB.ArcLength(0_s,min) / 1_m);
- histBlog(lineWithB.ArcLength(0_s, min) / 1_m);
- }
-
- int pdg = static_cast<int>(particles::GetPDG(p.GetPID()));
- if (abs(pdg) == 13) {
- histLmugeo(L);
- histELSrelmu(L * 1_m / lineWithB.ArcLength(0_s, min) - 1, p.GetEnergy() / 1_GeV);
- }
- if (abs(pdg) == 211 || abs(pdg) == 111) {
- histLpigeo(L);
- if (chargeNumber != 0)
- histELSrelpi(L * 1_m / lineWithB.ArcLength(0_s, min) - 1, p.GetEnergy() / 1_GeV);
- }
- if (abs(pdg) == 2212 || abs(pdg) == 2112) {
- histLpgeo(L);
- if (chargeNumber != 0)
- histELSrelp(L * 1_m / lineWithB.ArcLength(0_s, min) - 1, p.GetEnergy() / 1_GeV);
- }
- if (abs(pdg) == 11) {
- histLegeo(L);
- if (chargeNumber != 0)
- histELSrele(L * 1_m / lineWithB.ArcLength(0_s, min) - 1, p.GetEnergy() / 1_GeV);
- }
- if (abs(pdg) == 22)
- histLymag(L);
-
-
- return std::make_tuple(geometry::Trajectory<geometry::Line>(line, min),
- velocity.norm() * min, Steplimit, minIter->second);
- }
-
- template <typename Particle> // was Stack previously, and argument was
- // Stack::StackIterator
-
- // determine direction of the particle after adding magnetic field
- corsika::geometry::Line MagneticStep(
- Particle const& p, corsika::geometry::Line line, corsika::units::si::LengthType Steplength, bool i) {
- using namespace corsika::units::si;
-
- if (line.GetV0().norm() * 1_s == 0_m)
- return line;
-
- //charge of the particle
- int chargeNumber;
- if (corsika::particles::IsNucleus(p.GetPID())) {
- chargeNumber = p.GetNuclearZ();
- }
- else {
- chargeNumber = corsika::particles::GetChargeNumber(p.GetPID());
- }
-
- auto const* currentLogicalVolumeNode = p.GetNode();
- auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(p.GetPosition());
- auto k = chargeNumber * corsika::units::constants::c * 1_eV / (p.GetMomentum().norm() * 1_V);
- geometry::Vector<dimensionless_d> direction = line.GetV0().normalized();
- auto position = p.GetPosition();
-
- // First Movement
- // assuming magnetic field does not change during movement
- position = position + direction * Steplength / 2;
- // Change of direction by magnetic field
- direction = direction + direction.cross(magneticfield) * Steplength * k;
- // Second Movement
- position = position + direction * Steplength / 2;
-
- if (i == true) {
- //if(chargeNumber != 0) {
- L = direction.norm() * Steplength / 2_m + Steplength / 2_m;
- //}
- }
-
- geometry::Vector<LengthType::dimension_type> const distance = position - p.GetPosition();
- if (distance.norm() == 0_m) {
- return line;
- }
- geometry::Line newLine(p.GetPosition(), distance.normalized() * line.GetV0().norm());
- return newLine;
- }
-
- template <typename Particle> // was Stack previously, and argument was
- // Stack::StackIterator
-
- // determine direction of the particle after adding magnetic field
- auto MagneticStep(Particle const& p, corsika::units::si::LengthType distance) {
- using namespace corsika::units::si;
-
- if (p.GetMomentum().norm() == 0_GeV) {
- return std::make_tuple(p.GetPosition(), p.GetMomentum() / 1_GeV, double(0));
- }
-
- //charge of the particle
- int chargeNumber;
- if (corsika::particles::IsNucleus(p.GetPID())) {
- chargeNumber = p.GetNuclearZ();
- }
- else {
- chargeNumber = corsika::particles::GetChargeNumber(p.GetPID());
- }
-
- auto const* currentLogicalVolumeNode = p.GetNode();
- auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(p.GetPosition());
- auto k = chargeNumber * corsika::units::constants::c * 1_eV / (p.GetMomentum().norm() * 1_V);
- geometry::Vector<dimensionless_d> direction = p.GetMomentum().normalized();
- auto position = p.GetPosition();
-
- // determine steplength for the magnetic field
- // because Steplength should not be distance
- LengthType Steplength = distance;
- if (chargeNumber != 0) {
- auto c = direction.cross(magneticfield) * chargeNumber * corsika::units::constants::c * 1_eV /
- (p.GetMomentum().norm() * 1_V);
- Steplength = sqrt(2 / c.squaredNorm() * (sqrt(c.squaredNorm() * distance * distance + 1) -1));
- std::cout << "Steplength " << Steplength << std::endl;
- }
- if (Steplength == 0_m) {
- Steplength = distance;
- }
- //This formula hasnt been tested properly
-
- // First Movement
- // assuming magnetic field does not change during movement
- position = position + direction * Steplength / 2;
- // Change of direction by magnetic field
- direction = direction + direction.cross(magneticfield) * Steplength * k;
- // Second Movement
- position = position + direction * Steplength / 2;
- auto L2 = direction.norm() * Steplength / 2_m + Steplength / 2_m;
- return std::make_tuple(position, direction.normalized(), L2);
- }
- };
-
- } // namespace tracking_line
-
-} // namespace corsika::process
-
-#endif
diff --git a/Processes/TrackingLine/TrackingLine_interpolation.h b/Processes/TrackingLine/TrackingLine_interpolation.h
deleted file mode 100644
index b712b6b06..000000000
--- a/Processes/TrackingLine/TrackingLine_interpolation.h
+++ /dev/null
@@ -1,193 +0,0 @@
-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#ifndef _include_corsika_processes_TrackingLine_h_
-#define _include_corsika_processes_TrackingLine_h_
-
-#include <corsika/geometry/Line.h>
-#include <corsika/geometry/Plane.h>
-#include <corsika/geometry/Sphere.h>
-#include <corsika/geometry/Trajectory.h>
-#include <corsika/geometry/Vector.h>
-#include <corsika/particles/ParticleProperties.h>
-#include <corsika/units/PhysicalUnits.h>
-#include <corsika/utl/quartic.h>
-#include <cmath>
-#include <optional>
-#include <type_traits>
-#include <utility>
-
-namespace corsika::environment {
- template <typename IEnvironmentModel>
- class Environment;
- template <typename IEnvironmentModel>
- class VolumeTreeNode;
-} // namespace corsika::environment
-
-namespace corsika::process {
-
- namespace tracking_line {
-
- std::optional<std::pair<corsika::units::si::TimeType, corsika::units::si::TimeType>>
- TimeOfIntersection(geometry::Line const&, geometry::Sphere const&);
-
- corsika::units::si::TimeType TimeOfIntersection(geometry::Line const&,
- geometry::Plane const&);
-
- class TrackingLine {
-
- public:
- TrackingLine(){};
-
- template <typename Particle> // was Stack previously, and argument was
- // Stack::StackIterator
- auto GetTrack(Particle const& p) {
- using namespace corsika::units::si;
- using namespace corsika::geometry;
- geometry::Vector<SpeedType::dimension_type> velocity =
- p.GetMomentum() / p.GetEnergy() * corsika::units::constants::c;
-
- auto const currentPosition = p.GetPosition();
- std::cout << "TrackingLine pid: " << p.GetPID()
- << " , E = " << p.GetEnergy() / 1_GeV << " GeV" << std::endl;
- std::cout << "TrackingLine pos: "
- << currentPosition.GetCoordinates()
- // << " [" << p.GetNode()->GetModelProperties().GetName() << "]"
- << std::endl;
- std::cout << "TrackingLine E: " << p.GetEnergy() / 1_GeV << " GeV" << std::endl;
-
- // determine velocity after adding magnetic field
- auto const* currentLogicalVolumeNode = p.GetNode();
- int chargeNumber;
- if(corsika::particles::IsNucleus(p.GetPID())) {
- chargeNumber = p.GetNuclearZ();
- } else {
- chargeNumber = corsika::particles::GetChargeNumber(p.GetPID());
- }
- geometry::Vector<dimensionless_d> const directionBefore = velocity.normalized();
- auto magMaxLength = 1_m/0;
- auto directionAfter = directionBefore;
- if(chargeNumber != 0) {
- auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(currentPosition);
- std::cout << "TrackingLine B: " << magneticfield.GetComponents() / 1_uT << " uT " << std::endl;
- geometry::Vector<SpeedType::dimension_type> const velocityVerticalMag = velocity -
- velocity.parallelProjectionOnto(magneticfield);
- LengthType const gyroradius = p.GetEnergy() * velocityVerticalMag.GetNorm() * 1_V /
- (corsika::units::constants::cSquared * abs(chargeNumber) *
- magneticfield.GetNorm() * 1_eV);
- //steplength should consider more things than just gyroradius
- double maxAngle = 1e-5;
- LengthType const Steplength = 2 * sin(maxAngle * M_PI / 180) * gyroradius;
-
- std::cout << "Test: " << Steplength << " " << gyroradius << std::endl;
-
- // First Movement
- auto position = currentPosition + directionBefore * Steplength / 2;
- // Change of direction by magnetic field at position
- magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(position);
- directionAfter = directionBefore + directionBefore.cross(magneticfield) * chargeNumber *
- Steplength * corsika::units::constants::cSquared * 1_eV /
- (p.GetEnergy() * velocity.GetNorm() * 1_V);
- // Second Movement
- position = position + directionAfter * Steplength / 2;
- magMaxLength = (position - currentPosition).GetNorm();
- geometry::Vector<dimensionless_d> const direction = (position - currentPosition) /
- magMaxLength;
- velocity = direction * velocity.GetNorm();
- std::cout << "TrackingLine p: " << (direction * p.GetMomentum().GetNorm()).GetComponents() / 1_GeV
- << " GeV " << std::endl;
-
- } else {
- std::cout << "TrackingLine p: " << p.GetMomentum().GetComponents() / 1_GeV
- << " GeV " << std::endl;
- }
-
- std::cout << "TrackingLine v: " << velocity.GetComponents() << std::endl;
-
- geometry::Line line(currentPosition, velocity);
-
- //auto const* currentLogicalVolumeNode = p.GetNode();
- //~ auto const* currentNumericalVolumeNode =
- //~ fEnvironment.GetUniverse()->GetContainingNode(currentPosition);
- auto const numericallyInside =
- currentLogicalVolumeNode->GetVolume().Contains(currentPosition);
-
- std::cout << "numericallyInside = " << (numericallyInside ? "true" : "false");
-
- auto const& children = currentLogicalVolumeNode->GetChildNodes();
- auto const& excluded = currentLogicalVolumeNode->GetExcludedNodes();
-
- std::vector<std::pair<TimeType, decltype(p.GetNode())>> intersections;
-
- // for entering from outside
- auto addIfIntersects = [&](auto const& vtn) {
- auto const& volume = vtn.GetVolume();
- auto const& sphere = dynamic_cast<geometry::Sphere const&>(
- volume); // for the moment we are a bit bold here and assume
- // everything is a sphere, crashes with exception if not
-
- std::cout << "Mittelpunkt: " << sphere.GetCenter().GetCoordinates() << std::endl;
-
- if (auto opt = TimeOfIntersection(line, sphere); opt.has_value()) {
- auto const [t1, t2] = *opt;
- std::cout << "intersection times: " << t1 / 1_s << "; "
- << t2 / 1_s
- // << " " << vtn.GetModelProperties().GetName()
- << std::endl;
- if (t1.magnitude() > 0)
- intersections.emplace_back(t1, &vtn);
- else if (t2.magnitude() > 0)
- std::cout << "inside other volume" << std::endl;
- }
- };
-
- for (auto const& child : children) { addIfIntersects(*child); }
- for (auto const* ex : excluded) { addIfIntersects(*ex); }
-
- {
- auto const& sphere = dynamic_cast<geometry::Sphere const&>(
- currentLogicalVolumeNode->GetVolume());
- // for the moment we are a bit bold here and assume
- // everything is a sphere, crashes with exception if not
- [[maybe_unused]] auto const [t1, t2] = *TimeOfIntersection(line, sphere);
- [[maybe_unused]] auto dummy_t1 = t1;
- intersections.emplace_back(t2, currentLogicalVolumeNode->GetParent());
- }
-
- auto const minIter = std::min_element(
- intersections.cbegin(), intersections.cend(),
- [](auto const& a, auto const& b) { return a.first < b.first; });
-
- TimeType min;
-
- if (minIter == intersections.cend()) {
- min = 1_s; // todo: do sth. more reasonable as soon as tracking is able
- // to handle the numerics properly
- throw std::runtime_error("no intersection with anything!");
- } else {
- min = minIter->first;
- }
-
- std::cout << " t-intersect: "
- << min
- // << " " << minIter->second->GetModelProperties().GetName()
- << std::endl;
-
- return std::make_tuple(geometry::Trajectory<geometry::Line>(line, min),
- velocity.norm() * min, minIter->second, magMaxLength,
- directionBefore, directionAfter);
- }
- };
-
- } // namespace tracking_line
-
-} // namespace corsika::process
-
-#endif
diff --git a/Processes/TrackingLine/dump_bh.hpp b/Processes/TrackingLine/dump_bh.hpp
deleted file mode 100644
index 97f8460d7..000000000
--- a/Processes/TrackingLine/dump_bh.hpp
+++ /dev/null
@@ -1,186 +0,0 @@
-#pragma once
-
-#include <boost/histogram.hpp>
-#include <sstream>
-#include <string>
-#include <vector>
-
-/*
- To save one 1D boost::histogram in C++, include this file and run:
-
- dump_bh_1d("file_name.json", hist);
-
- Note that the file_name.json will be overwritten, so use a new file for each histogram! or use:
-
- std::ofstream file1("test_hist.json");
- dump_bh(file1, h1);
- file << ",\n";
- dump_bh(file1, h2);
- file << ",\n";
- dump_bh(file1, h3);
- file1.close();
-
-
- In python, just read the json and plot this histogram with matplotib
-
- for 1D histogram:
-
- import json
- import matplotlib.pyplot as plt
-
- h1=json.load(open("test_hist_direct.json"))
- print (h1)
- plt.hist(h1['bins'][:-1], bins=h1['bins'], weights=h1['data'])
- plt.show()
-
-
- for 2D histogram (needs some list/array gymnastics):
-
- import json
- import matplotlib.pyplot as plt
- import numpy as np
-
- h2=json.load(open("test_hist_2D_direct.json"))
- print (h2)
- xx,yy = np.meshgrid(h2['xbins'][:-1], h2['ybins'][:-1])
- plt.hist2d([i for s in xx for i in s], [i for s in yy for i in s],
- bins=[h2['xbins'],h2['ybins']], weights=h2['data'])
- plt.show()
-
-
- Do not forget to add axis titles, legend, etc.
- */
-
-
-template<typename T>
-void dump_bh_2d(std::ostream& os, T& h)
-{
- // determine number of axes (rank) and axes sizes
- const int rank = 2;
- std::vector<int> axis_size(rank);
- for (unsigned int i=0; i<rank; ++i)
- axis_size[i] = h.axis(i).size();
-
- // start dump json object
- os << "{\n";
- os << " \"rank\" : " << rank << ",\n";
- os << " \"size\" : [";
- bool first = true;
- for (auto s : axis_size) {
- os << (first?"":", ") << s;
- first = false;
- }
- os << "],\n";
-
- // bins-x
- os << " \"xbins\" : [";
- double lastx = 0;
- for (int i=0; i<h.axis(0).size(); ++i) {
- os << h.axis(0).bin(i).lower() << ", ";
- lastx = h.axis(0).bin(i).upper();
- }
- os << lastx << "],\n";
-
- // bins-y
- os << " \"ybins\" : [";
- double lasty = 0;
- for (int i=0; i<h.axis(1).size(); ++i) {
- os << h.axis(1).bin(i).lower() << ", ";
- lasty = h.axis(1).bin(i).upper();
- }
- os << lasty << "],\n";
-
- // binned data
- os << " \"data\" : [";
- first = true;
- for (auto x : boost::histogram::indexed(h)) {
- os << (first?"":", ") << *x;
- first = false;
- }
- os << " ]\n";
- os << "}\n";
-}
-
-
-
-
-template<typename T>
-void dump_bh_1d(std::ostream& os, T& h)
-{
- const int rank = 1;
- // determine number of axes (rank) and axes sizes
- std::vector<int> axis_size(rank);
- for (unsigned int i=0; i<rank; ++i)
- axis_size[i] = h.axis(i).size();
-
- // start dump json object
- os << "{\n";
- os << " \"rank\" : 1,\n";
- os << " \"size\" : [";
- bool first = true;
- for (auto s : axis_size) {
- os << (first?"":", ") << s << "],\n";
- first = false;
- }
-
- // underflow data
- os << " \"underflow\" : [";
- os << h.at(boost::histogram::axis::option::underflow);
- os << "],\n";
-
- // overflow data
- os << " \"overflow\" : [";
- os << h.at(boost::histogram::axis::option::overflow);
- os << "],\n";
-
- // bins
- os << " \"bins\" : [";
- double last = 0;
- for (auto x : boost::histogram::indexed(h)) {
- os << x.bin().lower() << ", ";
- last = x.bin().upper();
- }
- os << last << "],\n";
-
- // data of bins
- os << " \"data\" : [";
- first = true;
- for (auto x : boost::histogram::indexed(h)) {
- os << (first?"":", ") << *x;
- first = false;
- }
- os << " ]\n";
- os << "}\n";
-}
-
-
-
-
-template<typename T>
-void dump_bh(std::ostream& os, T& h)
-{
- // determine number of axes (rank) and axes sizes
- const int rank = h.rank();
- switch (rank) {
- case 1:
- dump_bh_1d(os, h);
- break;
- case 2:
- dump_bh_2d(os, h);
- break;
- default:
- {
- throw std::runtime_error("multi dim hist not implemented");
- }
- break;
- }
-}
-
-template<typename T>
-void dump_bh(const char* fname, T& h)
-{
- std::ofstream file(fname);
- dump_bh(file, h);
- file.close();
-}
-
diff --git a/Processes/TrackingLine/testTrackingLine.cc b/Processes/TrackingLine/testTrackingLine.cc
index 7634758e9..df6524185 100644
--- a/Processes/TrackingLine/testTrackingLine.cc
+++ b/Processes/TrackingLine/testTrackingLine.cc
@@ -6,7 +6,7 @@
* the license.
*/
-#include <corsika/process/tracking_line/TrackingLine.h>
+#include <corsika/process/tracking_line/Tracking.h>
#include <testTrackingLineStack.h> // test-build, and include file is obtained from CMAKE_CURRENT_SOURCE_DIR
#include <corsika/environment/Environment.h>
@@ -15,6 +15,7 @@
#include <corsika/geometry/Point.h>
#include <corsika/geometry/Sphere.h>
#include <corsika/geometry/Vector.h>
+#include <corsika/geometry/Intersections.hpp>
#include <corsika/setup/SetupTrajectory.h>
using corsika::setup::Trajectory;
@@ -30,74 +31,5 @@ using namespace corsika::geometry;
using namespace std;
using namespace corsika::units::si;
-
TEST_CASE("TrackingLine") {
- environment::Environment<TestMagneticField> env; // dummy environment
- auto const& cs = env.GetCoordinateSystem();
-
- tracking_line::TrackingLine tracking;
-
- SECTION("intersection with sphere") {
- Point const origin(cs, {0_m, 0_m, -5_m});
- Point const center(cs, {0_m, 0_m, 10_m});
- Sphere const sphere(center, 1_m);
- Vector<corsika::units::si::SpeedType::dimension_type> v(cs, 0_m / second,
- 0_m / second, 1_m / second);
- Line line(origin, v);
-
- setup::Trajectory traj(line, 12345_s);
-
- auto const opt =
- tracking_line::TimeOfIntersection(traj, Sphere(Point(cs, {0_m, 0_m, 10_m}), 1_m));
- REQUIRE(opt.has_value());
-
- auto [t1, t2] = opt.value();
- REQUIRE(t1 / 14_s == Approx(1));
- REQUIRE(t2 / 16_s == Approx(1));
-
- auto const optNoIntersection =
- tracking_line::TimeOfIntersection(traj, Sphere(Point(cs, {5_m, 0_m, 10_m}), 1_m));
- REQUIRE_FALSE(optNoIntersection.has_value());
- }
-
- SECTION("maximally possible propagation") {
- auto& universe = *(env.GetUniverse());
-
- auto const radius = 20_m;
-
- auto theMedium = environment::Environment<TestMagneticField>::CreateNode<Sphere>(
- Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m}, radius);
- auto const* theMediumPtr = theMedium.get();
- universe.AddChild(std::move(theMedium));
-
- TestTrackingLineStack stack;
- stack.AddParticle(
- std::tuple<particles::Code, units::si::HEPEnergyType,
- corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
- particles::Code::MuPlus,
- 1_GeV,
- {cs, {0_GeV, 0_GeV, 1_GeV}},
- {cs, {0_m, 0_m, 0_km}},
- 0_ns});
- auto p = stack.GetNextParticle();
- p.SetNode(theMediumPtr);
-
- Point const origin(cs, {0_m, 0_m, 0_m});
- Vector<corsika::units::si::SpeedType::dimension_type> v(cs, 0_m / second,
- 0_m / second, 1_m / second);
- Line line(origin, v);
-
- const auto [stepWithoutB, stepWithB, geomMaxLength, magMaxLength, nextVol] = tracking.GetTrack(p);
- //auto const [traj, geomMaxLength, nextVol, magMaxLength, beforeDirection,
- // afterDirection] = tracking.GetTrack(p);
- [[maybe_unused]] auto& dummy_1 = stepWithB;
- [[maybe_unused]] auto& dummy_2 = magMaxLength;
- [[maybe_unused]] auto& dummy_geomMaxLength = geomMaxLength;
- [[maybe_unused]] auto& dummy_nextVol = nextVol;
-
- REQUIRE((stepWithoutB.GetPosition(1.) - Point(cs, 0_m, 0_m, radius))
- .GetComponents(cs)
- .norm()
- .magnitude() == Approx(0).margin(1e-4));
- }
}
diff --git a/Processes/TrackingLine/testTrackingLineStack.h b/Processes/TrackingLine/testTrackingLineStack.h
index 16b07e0cb..5c714048a 100644
--- a/Processes/TrackingLine/testTrackingLineStack.h
+++ b/Processes/TrackingLine/testTrackingLineStack.h
@@ -21,21 +21,8 @@
#include <corsika/units/PhysicalUnits.h>
-
-class TestMagneticField {
- using MagneticFieldVector =
- corsika::geometry::Vector<corsika::units::si::magnetic_flux_density_d>;
-
-public:
- MagneticFieldVector GetMagneticField(corsika::geometry::Point const& p) const {
- using namespace corsika::units::si;
- return MagneticFieldVector(p.GetCoordinateSystem(), 0_T, 0_T, 1_T);
- }
-};
-
-
using TestEnvironmentType =
- corsika::environment::Environment<TestMagneticField>;
+ corsika::environment::Environment<corsika::environment::Empty>;
template <typename T>
using SetupGeometryDataInterface =
diff --git a/Setup/SetupEnvironment.h b/Setup/SetupEnvironment.h
index 9c893beb2..87f0c8e30 100644
--- a/Setup/SetupEnvironment.h
+++ b/Setup/SetupEnvironment.h
@@ -40,7 +40,9 @@ namespace corsika::setup {
*/
namespace corsika::setup::testing {
- inline auto setupEnvironment(particles::Code vTargetCode) {
+ inline auto setupEnvironment(particles::Code vTargetCode,
+ const corsika::units::si::MagneticFluxType BfieldZ =
+ corsika::units::si::MagneticFluxType::zero()) {
using namespace corsika::units::si;
using namespace corsika;
@@ -53,8 +55,7 @@ namespace corsika::setup::testing {
* our world is a sphere at 0,0,0 with R=infty
*/
auto world = setup::Environment::CreateNode<geometry::Sphere>(
- geometry::Point{cs, 0_m, 0_m, 0_m},
- 1_km * std::numeric_limits<double>::infinity());
+ geometry::Point{cs, 0_m, 0_m, 0_m}, 100_km);
/**
* construct suited environment medium model:
@@ -64,12 +65,12 @@ namespace corsika::setup::testing {
environment::HomogeneousMedium<setup::EnvironmentInterface>>>;
world->SetModelProperties<MyHomogeneousModel>(
- environment::Medium::AirDry1Atm, geometry::Vector(cs, 0_T, 0_T, 1_T),
+ environment::Medium::AirDry1Atm, geometry::Vector(cs, 0_T, 0_T, BfieldZ),
1_kg / (1_m * 1_m * 1_m),
environment::NuclearComposition(std::vector<particles::Code>{vTargetCode},
std::vector<float>{1.}));
- auto const* nodePtr = world.get();
+ auto* nodePtr = world.get();
universe.AddChild(std::move(world));
return std::make_tuple(std::move(env), &cs, nodePtr);
diff --git a/Setup/SetupTrajectory.h b/Setup/SetupTrajectory.h
index 53fe4d04f..a441bfcab 100644
--- a/Setup/SetupTrajectory.h
+++ b/Setup/SetupTrajectory.h
@@ -19,40 +19,6 @@
namespace corsika::setup {
/// definition of Trajectory base class, to be used in tracking and cascades
- typedef corsika::geometry::Trajectory<corsika::geometry::Line> Trajectory;
+ typedef corsika::geometry::LineTrajectory Trajectory;
- /*
- typedef std::variant<std::monostate, corsika::geometry::Trajectory<Line>,
- corsika::geometry::Trajectory<Helix>>
- Trajectory;
-
- /// helper visitor to modify Particle by moving along Trajectory
- template <typename Particle>
- class ParticleUpdate {
-
- Particle& fP;
-
- public:
- ParticleUpdate(Particle& p)
- : fP(p) {}
- void operator()(std::monostate const&) {}
-
- template <typename T>
- void operator()(T const& trajectory) {
- fP.SetPosition(trajectory.GetPosition(1));
- }
- };
-
- /// helper visitor to modify Particle by moving along Trajectory
- class GetDuration {
- public:
- corsika::units::si::TimeType operator()(std::monostate const&) {
- return 0 * corsika::units::si::second;
- }
- template <typename T>
- corsika::units::si::TimeType operator()(T const& trajectory) {
- return trajectory.GetDuration();
- }
- };
- */
} // namespace corsika::setup
diff --git a/Stack/History/HistoryObservationPlane.cpp b/Stack/History/HistoryObservationPlane.cpp
index b1425f18b..eb5092d1f 100644
--- a/Stack/History/HistoryObservationPlane.cpp
+++ b/Stack/History/HistoryObservationPlane.cpp
@@ -28,12 +28,12 @@ HistoryObservationPlane::HistoryObservationPlane(setup::Stack const& stack,
corsika::process::EProcessReturn HistoryObservationPlane::DoContinuous(
setup::Stack::ParticleType const& particle, setup::Trajectory const& trajectory) {
TimeType const timeOfIntersection =
- (plane_.GetCenter() - trajectory.GetR0()).dot(plane_.GetNormal()) /
- trajectory.GetV0().dot(plane_.GetNormal());
+ (plane_.GetCenter() - trajectory.GetLine().GetR0()).dot(plane_.GetNormal()) /
+ trajectory.GetLine().GetV0().dot(plane_.GetNormal());
if (timeOfIntersection < TimeType::zero()) { return process::EProcessReturn::eOk; }
- if (plane_.IsAbove(trajectory.GetR0()) == plane_.IsAbove(trajectory.GetPosition(1))) {
+ if (plane_.IsAbove(trajectory.GetLine().GetR0()) == plane_.IsAbove(trajectory.GetPosition(1))) {
return process::EProcessReturn::eOk;
}
@@ -53,15 +53,15 @@ corsika::process::EProcessReturn HistoryObservationPlane::DoContinuous(
LengthType HistoryObservationPlane::MaxStepLength(setup::Stack::ParticleType const&,
setup::Trajectory const& trajectory) {
TimeType const timeOfIntersection =
- (plane_.GetCenter() - trajectory.GetR0()).dot(plane_.GetNormal()) /
- trajectory.GetV0().dot(plane_.GetNormal());
+ (plane_.GetCenter() - trajectory.GetLine().GetR0()).dot(plane_.GetNormal()) /
+ trajectory.GetLine().GetV0().dot(plane_.GetNormal());
if (timeOfIntersection < TimeType::zero()) {
return std::numeric_limits<double>::infinity() * 1_m;
}
- auto const pointOfIntersection = trajectory.GetPosition(timeOfIntersection);
- return (trajectory.GetR0() - pointOfIntersection).norm() * 1.0001;
+ auto const pointOfIntersection = trajectory.GetLine().GetPosition(timeOfIntersection);
+ return (trajectory.GetLine().GetR0() - pointOfIntersection).norm() * 1.0001;
}
void HistoryObservationPlane::fillHistoryHistogram(
diff --git a/do-copyright.py b/do-copyright.py
index 716b2c4bc..6b235afce 100755
--- a/do-copyright.py
+++ b/do-copyright.py
@@ -21,10 +21,10 @@ Debug settings are 0: nothing, 1: checking, 2: filesystem
"""
Debug = 0
-excludeDirs = ["ThirdParty", "git", "build", "install", "PROPOSAL"]
+excludeDirs = ["modules", "git", "build", "install", "externals"]
excludeFiles = ['PhysicalConstants.h','CorsikaFenvOSX.cc', 'sgn.h', 'quartic.h']
-extensions = [".cc", ".h", ".test"]
+extensions = [".cpp", ".hpp"]
"""
justCheck: T: only checking, F: also changing files
@@ -192,9 +192,10 @@ def next_file(dir_name, files, justCheck, forYear, updateMessage):
excludes if wished, process otherwise
"""
for check in excludeDirs :
+ print (dir_name)
if check in dir_name:
- if Debug>1:
- print ("exclude-dir: " + check)
+ # if Debug>1:
+ print ("exclude-dir: " + check)
return True
for check in files :
if (os.path.isdir(check)):
--
GitLab