diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 464ce2402377ddacbf94a9147aa958ebbd9bdf3d..9ebe0555b9b6f4745c25e3e1bf09a3c1723a05eb 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -170,7 +170,7 @@ build-clang-8:
script:
- cd build
- set -o pipefail
- - ctest -j4 -VV | gzip -v -9 > test.log.gz
+ - ctest -j4
rules:
- if: $CI_MERGE_REQUEST_ID
when: manual
@@ -185,8 +185,6 @@ build-clang-8:
reports:
junit:
- ${CI_PROJECT_DIR}/build/test_outputs/junit*.xml
- paths:
- - ${CI_PROJECT_DIR}/build/test.log.gz
cache:
paths:
- ${CI_PROJECT_DIR}/build/
@@ -227,7 +225,7 @@ test-clang-8:
- cd build
- cmake --build . -- -j4
- set -o pipefail
- - ctest -j4 -VV | gzip -v -9 > test.log.gz
+ - ctest -j4
rules:
- if: '$CI_MERGE_REQUEST_ID && $CI_MERGE_REQUEST_TITLE =~ /^Draft:/'
allow_failure: false
@@ -241,8 +239,6 @@ test-clang-8:
reports:
junit:
- ${CI_PROJECT_DIR}/build/test_outputs/junit*.xml
- paths:
- - ${CI_PROJECT_DIR}/build/test.log.gz
cache:
paths:
- ${CI_PROJECT_DIR}/build/
@@ -338,7 +334,7 @@ example-clang-8:
- cd build
- cmake --build . -- -j4
- set -o pipefail
- - ctest -j4 -VV | gzip -v -9 > test.log.gz
+ - ctest -j4
- make -j4 run_examples | gzip -v -9 > examples.log.gz
rules:
- if: '$CI_MERGE_REQUEST_ID && $CI_MERGE_REQUEST_TITLE =~ /^Draft:/'
@@ -357,7 +353,6 @@ example-clang-8:
- ${CI_PROJECT_DIR}/build/test_outputs/junit*.xml
paths:
- ${CI_PROJECT_DIR}/build/examples.log.gz
- - ${CI_PROJECT_DIR}/build/test.log.gz
cache:
paths:
- ${CI_PROJECT_DIR}/build/
@@ -450,7 +445,7 @@ install-clang-8:
- cmake .. -DCMAKE_BUILD_TYPE=Release
- cmake --build . -- -j4
- set -o pipefail
- - ctest -j4 -VV | gzip -v -9 > test.log.gz
+ - ctest -j4
- make -j4 run_examples | gzip -v -9 > examples.log.gz
rules:
- if: '$CI_MERGE_REQUEST_LABELS =~ /Ready for code review/' # run on merge requests, if label 'Ready for code review' is set
@@ -475,7 +470,6 @@ install-clang-8:
junit:
- ${CI_PROJECT_DIR}/build/test_outputs/junit*.xml
paths:
- - ${CI_PROJECT_DIR}/build/test.log.gz
- ${CI_PROJECT_DIR}/build/examples.log.gz
# release for gcc
@@ -513,7 +507,7 @@ coverage:
- cd build
- cmake .. -DCMAKE_BUILD_TYPE=Coverage
- cmake --build . -- -j4
- - ctest -j4 -VV | gzip -v -9 > test.log.gz
+ - ctest -j4
- cmake --build . --target coverage
- tar czf coverage-report.tar.gz coverage-report
coverage: '/^.*functions\.+:\s(.*\%)\s/'
@@ -530,7 +524,6 @@ coverage:
when: always
expire_in: 1 year
paths:
- - ${CI_PROJECT_DIR}/build/test.log.gz
- ${CI_PROJECT_DIR}/build/coverage-report.tar.gz
cache:
paths:
diff --git a/CMakeLists.txt b/CMakeLists.txt
index b126b3eed5da47088f15ed56d50078c6aa9dbbe2..e24a55b9592c0917db26acdab4fbdaf89a27d469 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/Data b/Data
index 150c0b40b0fe4c44117d064f8184de864d4a958a..4d1785e74a7bf6769ee016509b357e63d7117438 160000
--- a/Data
+++ b/Data
@@ -1 +1 @@
-Subproject commit 150c0b40b0fe4c44117d064f8184de864d4a958a
+Subproject commit 4d1785e74a7bf6769ee016509b357e63d7117438
diff --git a/Documentation/Examples/boundary_example.cc b/Documentation/Examples/boundary_example.cc
index c8fa31e6edf2c0457e5fe6617897656bac5bb5f4..ed9d7bd43cc24c03df3642315cd484afde1fb99b 100644
--- a/Documentation/Examples/boundary_example.cc
+++ b/Documentation/Examples/boundary_example.cc
@@ -8,7 +8,6 @@
#include <corsika/cascade/Cascade.h>
#include <corsika/process/ProcessSequence.h>
-#include <corsika/process/tracking_line/TrackingLine.h>
#include <corsika/setup/SetupEnvironment.h>
#include <corsika/setup/SetupStack.h>
@@ -85,7 +84,7 @@ private:
//
int main() {
- logging::SetLevel(logging::level::info);
+ logging::SetLevel(logging::level::trace);
C8LOG_INFO("boundary_example");
@@ -102,7 +101,7 @@ int main() {
// create "world" as infinite sphere filled with protons
auto world = EnvType::CreateNode<Sphere>(
- Point{rootCS, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
+ Point{rootCS, 0_m, 0_m, 0_m}, 100_km);
using MyHomogeneousModel =
environment::MediumPropertyModel<environment::UniformMagneticField<
@@ -123,7 +122,7 @@ int main() {
universe.AddChild(std::move(world));
// setup processes, decays and interactions
- tracking_line::TrackingLine tracking;
+ setup::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 f3cd9c19aa48d5eb02c7276d10d07752f2a2cfbc..8b9b5b8ad0b0c50830b01ec21fd7dc6163c91d20 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -10,7 +10,6 @@
#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/setup/SetupEnvironment.h>
#include <corsika/setup/SetupStack.h>
@@ -58,6 +57,8 @@ int main() {
logging::SetLevel(logging::level::info);
+ C8LOG_INFO("vertical_EAS");
+
std::cout << "cascade_example" << std::endl;
const LengthType height_atmosphere = 112.8_km;
@@ -73,7 +74,7 @@ int main() {
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
auto world = setup::Environment::CreateNode<Sphere>(
- Point{rootCS, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
+ Point{rootCS, 0_m, 0_m, 0_m}, 150_km);
using MyHomogeneousModel =
environment::MediumPropertyModel<environment::UniformMagneticField<
@@ -111,7 +112,7 @@ int main() {
rootCS, 0_m, 0_m,
height_atmosphere); // this is the CORSIKA 7 start of atmosphere/universe
- ShowerAxis const showerAxis{injectionPos, Vector{rootCS, 0_m, 0_m, -5000_km}, env};
+ ShowerAxis const showerAxis{injectionPos, Vector{rootCS, 0_m, 0_m, -100_km}, env};
{
auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
@@ -135,7 +136,7 @@ int main() {
}
// setup processes, decays and interactions
- tracking_line::TrackingLine tracking;
+ setup::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 2c5e1489e4fe9baa98d6e879e3fe1fee0b70ee79..b36239219320f05e006fe078404f5e0d51eef117 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/energy_loss/EnergyLoss.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
@@ -18,6 +18,7 @@
#include <corsika/environment/Environment.h>
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/NuclearComposition.h>
+#include <corsika/environment/ShowerAxis.h>
#include <corsika/geometry/Sphere.h>
@@ -74,20 +75,20 @@ int main() {
auto& universe = *(env.GetUniverse());
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
- auto theMedium = EnvType::CreateNode<Sphere>(
- Point{rootCS, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
+ auto world = EnvType::CreateNode<Sphere>(
+ Point{rootCS, 0_m, 0_m, 0_m}, 150_km);
using MyHomogeneousModel =
environment::MediumPropertyModel<environment::UniformMagneticField<
environment::HomogeneousMedium<setup::EnvironmentInterface>>>;
- theMedium->SetModelProperties<MyHomogeneousModel>(
+ world->SetModelProperties<MyHomogeneousModel>(
environment::Medium::AirDry1Atm, geometry::Vector(rootCS, 0_T, 0_T, 1_T),
1_kg / (1_m * 1_m * 1_m),
NuclearComposition(std::vector<particles::Code>{particles::Code::Hydrogen},
std::vector<float>{(float)1.}));
- universe.AddChild(std::move(theMedium));
+ universe.AddChild(std::move(world));
// setup particle stack, and add primary particle
setup::Stack stack;
@@ -98,6 +99,7 @@ int main() {
double theta = 0.;
double phi = 0.;
+ Point injectionPos(rootCS, 0_m, 0_m, 0_m);
{
auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
return sqrt(Elab * Elab - m * m);
@@ -113,15 +115,14 @@ int main() {
cout << "input particle: " << beamCode << endl;
cout << "input angles: theta=" << theta << " phi=" << phi << endl;
cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl;
- Point pos(rootCS, 0_m, 0_m, 0_m);
stack.AddParticle(
std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
- beamCode, E0, plab, pos, 0_ns});
+ beamCode, E0, plab, injectionPos, 0_ns});
}
// setup processes, decays and interactions
- tracking_line::TrackingLine tracking;
+ setup::Tracking tracking;
stack_inspector::StackInspector<setup::Stack> stackInspect(1, true, E0);
random::RNGManager::GetInstance().RegisterRandomStream("sibyll");
@@ -131,20 +132,19 @@ int main() {
// process::sibyll::NuclearInteraction sibyllNuc(env, sibyll);
// process::sibyll::Decay decay;
process::pythia::Decay decay;
- process::particle_cut::ParticleCut cut(20_GeV, true, true);
+ process::particle_cut::ParticleCut cut(60_GeV, true, true);
// random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel");
// process::HadronicElasticModel::HadronicElasticInteraction
// hadronicElastic(env);
process::track_writer::TrackWriter trackWriter("tracks.dat");
+ ShowerAxis const showerAxis{injectionPos, Vector{rootCS, 0_m, 0_m, -100_km}, env};
+ process::energy_loss::EnergyLoss eLoss{showerAxis, cut.GetECut()};
// assemble all processes into an ordered process list
// auto sequence = sibyll << decay << hadronicElastic << cut << trackWriter;
- auto sequence = process::sequence(pythia, decay, cut, trackWriter, stackInspect);
-
- // cout << "decltype(sequence)=" << type_id_with_cvr<decltype(sequence)>().pretty_name()
- // << "\n";
+ auto sequence = process::sequence(pythia, decay, eLoss, cut, trackWriter, stackInspect);
// define air shower object, run simulation
cascade::Cascade EAS(env, tracking, sequence, stack);
diff --git a/Documentation/Examples/em_shower.cc b/Documentation/Examples/em_shower.cc
index 594a95845c9a2e8d8640aeda4c7847ed591d9420..07c4cfb7b8dac86bb1dc389fc1c72ab896a00c75 100644
--- a/Documentation/Examples/em_shower.cc
+++ b/Documentation/Examples/em_shower.cc
@@ -21,7 +21,6 @@
#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/random/RNGManager.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
@@ -72,13 +71,13 @@ int main(int argc, char** argv) {
// setup environment, geometry
using EnvType = setup::Environment;
EnvType env;
- const CoordinateSystem& rootCS = env.GetCoordinateSystem();
+ const CoordinateSystem& rootCS = env.GetCoordinateSystem();
Point const center{rootCS, 0_m, 0_m, 0_m};
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});
+ 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
@@ -151,13 +150,13 @@ int main(int argc, char** argv) {
process::longitudinal_profile::LongitudinalProfile longprof{showerAxis};
Plane const obsPlane(showerCore, Vector<dimensionless_d>(rootCS, {0., 0., 1.}));
- process::observation_plane::ObservationPlane observationLevel(obsPlane,
- "particles.dat");
+ process::observation_plane::ObservationPlane observationLevel(
+ 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;
+ setup::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 c6647209bb13d4843020c2bce4d0354486f97753..2216e4b33b325b09bf4c209d154bf79ceda2bb66 100644
--- a/Documentation/Examples/hybrid_MC.cc
+++ b/Documentation/Examples/hybrid_MC.cc
@@ -34,7 +34,6 @@
#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/urqmd/UrQMD.h>
#include <corsika/random/RNGManager.h>
#include <corsika/setup/SetupStack.h>
@@ -102,7 +101,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
@@ -216,8 +215,8 @@ int main(int argc, char** argv) {
process::longitudinal_profile::LongitudinalProfile longprof{showerAxis};
Plane const obsPlane(showerCore, Vector<dimensionless_d>(rootCS, {0., 0., 1.}));
- process::observation_plane::ObservationPlane observationLevel(obsPlane,
- "particles.dat");
+ process::observation_plane::ObservationPlane observationLevel(
+ obsPlane, Vector<dimensionless_d>(rootCS, {1., 0., 0.}), "particles.dat");
process::UrQMD::UrQMD urqmd;
process::interaction_counter::InteractionCounter urqmdCounted{urqmd};
@@ -242,7 +241,7 @@ int main(int argc, char** argv) {
eLoss, cut, conex, longprof, observationLevel);
// define air shower object, run simulation
- tracking_line::TrackingLine tracking;
+ setup::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 c6f2dbe6111e7141f2880e320c83a4bdead4a030..6fbed4a73ecfd8e99db28d11754037e094033dc7 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -15,9 +15,13 @@
#include <corsika/cascade/Cascade.h>
#include <corsika/environment/Environment.h>
#include <corsika/environment/FlatExponential.h>
+#include <corsika/environment/HomogeneousMedium.h>
+#include <corsika/environment/IMagneticFieldModel.h>
#include <corsika/environment/LayeredSphericalAtmosphereBuilder.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/environment/ShowerAxis.h>
+#include <corsika/environment/SlidingPlanarExponential.h>
+#include <corsika/environment/UniformMagneticField.h>
#include <corsika/geometry/Plane.h>
#include <corsika/geometry/Sphere.h>
#include <corsika/logging/Logging.h>
@@ -29,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/stack_inspector/StackInspector.h>
#include <corsika/process/sibyll/Interaction.h>
#include <corsika/process/sibyll/NuclearInteraction.h>
-#include <corsika/process/tracking_line/TrackingLine.h>
#include <corsika/process/urqmd/UrQMD.h>
#include <corsika/random/RNGManager.h>
#include <corsika/setup/SetupStack.h>
@@ -105,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
@@ -148,6 +153,7 @@ int main(int argc, char** argv) {
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 +
@@ -159,8 +165,20 @@ int main(int argc, char** argv) {
stack.AddParticle(std::make_tuple(beamCode, E0, plab, injectionPos, 0_ns, A, Z));
} else {
- stack.AddParticle(
- std::make_tuple(particles::Code::Proton, E0, plab, injectionPos, 0_ns));
+ if (Z == 1) {
+ stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point,
+ units::si::TimeType>{particles::Code::Proton, E0, plab,
+ injectionPos, 0_ns});
+ } else if (Z == 0) {
+ stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point,
+ units::si::TimeType>{particles::Code::Neutron, E0,
+ plab, injectionPos, 0_ns});
+ } else {
+ std::cerr << "illegal parameters" << std::endl;
+ return EXIT_FAILURE;
+ }
}
// we make the axis much longer than the inj-core distance since the
@@ -173,6 +191,11 @@ int main(int argc, char** argv) {
// setup processes, decays and interactions
+ 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::sibyll::Interaction sibyll;
process::interaction_counter::InteractionCounter sibyllCounted(sibyll);
@@ -204,18 +227,14 @@ 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);
+ process::track_writer::TrackWriter trackWriter("tracks.dat");
process::longitudinal_profile::LongitudinalProfile longprof{showerAxis};
Plane const obsPlane(showerCore, Vector<dimensionless_d>(rootCS, {0., 0., 1.}));
- process::observation_plane::ObservationPlane observationLevel(obsPlane,
- "particles.dat");
+ process::observation_plane::ObservationPlane observationLevel(
+ obsPlane, Vector<dimensionless_d>(rootCS, {1., 0., 0.}), "particles.dat");
process::UrQMD::UrQMD urqmd;
process::interaction_counter::InteractionCounter urqmdCounted{urqmd};
@@ -236,12 +255,13 @@ int main(int argc, char** argv) {
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,
- proposalCounted, em_continuous, cut, observationLevel, longprof);
+ stack_inspector::StackInspector<setup::Stack> stackInspect(1000, false, E0);
+ auto sequence = process::sequence(stackInspect, hadronSequence, reset_particle_mass,
+ decaySequence, proposalCounted, em_continuous, cut,
+ trackWriter, observationLevel, longprof);
// define air shower object, run simulation
- tracking_line::TrackingLine tracking;
+ setup::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 97e2689bc249b79f33574149e4f9e2eb03c2e4eb..bf8c16d961c4cfdeb28827851dd8e20e8040c17e 100644
--- a/Environment/BaseExponential.h
+++ b/Environment/BaseExponential.h
@@ -10,9 +10,9 @@
#include <corsika/geometry/Line.h>
#include <corsika/geometry/Point.h>
-#include <corsika/geometry/Trajectory.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/units/PhysicalUnits.h>
+#include <corsika/setup/SetupTrajectory.h>
#include <limits>
@@ -47,12 +47,12 @@ namespace corsika::environment {
*/
// clang-format on
units::si::GrammageType IntegratedGrammage(
- geometry::Trajectory<geometry::Line> const& vLine, units::si::LengthType vL,
+ setup::Trajectory 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,10 @@ namespace corsika::environment {
*/
// clang-format on
units::si::LengthType ArclengthFromGrammage(
- geometry::Trajectory<geometry::Line> const& vLine,
- units::si::GrammageType vGrammage,
+ setup::Trajectory 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/CMakeLists.txt b/Environment/CMakeLists.txt
index c02aa95868b68b5c0304e1b76287b084a6b0debe..3231ab4a8a643b682a50723e9b64dc3837e1c761 100644
--- a/Environment/CMakeLists.txt
+++ b/Environment/CMakeLists.txt
@@ -19,6 +19,7 @@ set (
set (
ENVIRONMENT_HEADERS
VolumeTreeNode.h
+ IEmpty.hpp
IMediumModel.h
NuclearComposition.h
HomogeneousMedium.h
@@ -35,6 +36,7 @@ set (
ShowerAxis.h
IMagneticFieldModel.h
UniformMagneticField.h
+ NoMagneticField.h
IRefractiveIndexModel.h
UniformRefractiveIndex.h
IMediumPropertyModel.h
diff --git a/Environment/DensityFunction.h b/Environment/DensityFunction.h
index 5f61e26d7daf2502e47d8f6bc596d9cd484e5295..3ac33c3f21ebc0a8fbe792d9641b3a37c4dec7c4 100644
--- a/Environment/DensityFunction.h
+++ b/Environment/DensityFunction.h
@@ -11,7 +11,6 @@
#include <corsika/environment/LinearApproximationIntegrator.h>
#include <corsika/geometry/Line.h>
#include <corsika/geometry/Point.h>
-#include <corsika/geometry/Trajectory.h>
namespace corsika::environment {
diff --git a/Environment/Environment.h b/Environment/Environment.h
index 774b7492761469affbfc2d6bf326cac5501dbbf2..1f807de45bc698718c74f93cbc7021c8329936b2 100644
--- a/Environment/Environment.h
+++ b/Environment/Environment.h
@@ -8,7 +8,6 @@
#pragma once
-#include <corsika/environment/IMediumModel.h>
#include <corsika/environment/VolumeTreeNode.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/RootCoordinateSystem.h>
diff --git a/Environment/FlatExponential.h b/Environment/FlatExponential.h
index e1b0cd5b67c7976aee19f3642c02ea2c009b6ba6..8aeafbce4a776c1ae1dd4b335c3578962a21e215 100644
--- a/Environment/FlatExponential.h
+++ b/Environment/FlatExponential.h
@@ -12,10 +12,11 @@
#include <corsika/environment/NuclearComposition.h>
#include <corsika/geometry/Line.h>
#include <corsika/geometry/Point.h>
-#include <corsika/geometry/Trajectory.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/units/PhysicalUnits.h>
+#include <corsika/setup/SetupTrajectory.h>
+
namespace corsika::environment {
//clang-format off
@@ -51,14 +52,13 @@ namespace corsika::environment {
NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
- units::si::GrammageType IntegratedGrammage(
- geometry::Trajectory<geometry::Line> const& vLine,
- units::si::LengthType vTo) const override {
+ units::si::GrammageType IntegratedGrammage(setup::Trajectory const& vLine,
+ units::si::LengthType vTo) const override {
return Base::IntegratedGrammage(vLine, vTo, fAxis);
}
units::si::LengthType ArclengthFromGrammage(
- geometry::Trajectory<geometry::Line> const& vLine,
+ setup::Trajectory 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 49f4cb80d55f1fa229e0d324fa546db551a5f476..f39cc70610c7e9fdca38172e51f1d9c8179f494b 100644
--- a/Environment/HomogeneousMedium.h
+++ b/Environment/HomogeneousMedium.h
@@ -16,6 +16,8 @@
#include <corsika/random/RNGManager.h>
#include <corsika/units/PhysicalUnits.h>
+#include <corsika/setup/SetupTrajectory.h>
+
#include <cassert>
/**
@@ -42,14 +44,14 @@ namespace corsika::environment {
NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
corsika::units::si::GrammageType IntegratedGrammage(
- corsika::geometry::Trajectory<corsika::geometry::Line> const&,
+ corsika::setup::Trajectory 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::setup::Trajectory const&,
corsika::units::si::GrammageType pGrammage) const override {
return pGrammage / fDensity;
}
diff --git a/Environment/IEmpty.hpp b/Environment/IEmpty.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..6d45cd03c271868118d02d744fb4d6091e31cedf
--- /dev/null
+++ b/Environment/IEmpty.hpp
@@ -0,0 +1,47 @@
+/*
+ * (c) Copyright 2020 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 <corsika/geometry/Trajectory.h>
+
+namespace corsika::environment {
+
+ /**
+ * \class IEmpty
+ *
+ * intended for usage as default template argument for environments
+ * with no properties. For now, the ArclengthFromGrammage is
+ * mandatory, since it is used even in the most simple Cascade code.
+ *
+ * - IEmpty is the interface definition.
+ * - Empty<IEmpty> is a possible model implementation
+ *
+ **/
+
+ class IEmpty {
+ public:
+ virtual corsika::units::si::LengthType ArclengthFromGrammage(
+ corsika::geometry::LineTrajectory const&,
+ corsika::units::si::GrammageType) const = 0;
+
+ virtual ~IEmpty() {}
+ };
+
+ template <typename TModel = IEmpty>
+ class Empty : public TModel {
+ public:
+ corsika::units::si::LengthType ArclengthFromGrammage(
+ corsika::geometry::LineTrajectory const&,
+ corsika::units::si::GrammageType) const {
+ return 0. * corsika::units::si::meter;
+ }
+ };
+
+} // namespace corsika::environment
diff --git a/Environment/IMediumModel.h b/Environment/IMediumModel.h
index da6c9ca47a41ad14d58b1cc991a473ee970520bd..7a3b951711b8fe19f928eed3c5cf3d21a17e73d6 100644
--- a/Environment/IMediumModel.h
+++ b/Environment/IMediumModel.h
@@ -11,9 +11,10 @@
#include <corsika/environment/NuclearComposition.h>
#include <corsika/geometry/Line.h>
#include <corsika/geometry/Point.h>
-#include <corsika/geometry/Trajectory.h>
#include <corsika/units/PhysicalUnits.h>
+#include <corsika/setup/SetupTrajectory.h>
+
namespace corsika::environment {
class IMediumModel {
@@ -26,11 +27,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::setup::Trajectory const&,
corsika::units::si::LengthType) const = 0;
virtual corsika::units::si::LengthType ArclengthFromGrammage(
- corsika::geometry::Trajectory<corsika::geometry::Line> const&,
+ corsika::setup::Trajectory const&,
corsika::units::si::GrammageType) const = 0;
virtual NuclearComposition const& GetNuclearComposition() const = 0;
diff --git a/Environment/InhomogeneousMedium.h b/Environment/InhomogeneousMedium.h
index 60f5919235c690345c3a9c7f57a2b44af59cc8e9..76378d5d5202527b54d32811ec688f00807b0382 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::setup::Trajectory 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::setup::Trajectory 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 1224181f3a0e2d5fcb426eb7a70b6dd1f6a24026..3a6c37ccb6c7b92227093322a3024a864907fe0e 100644
--- a/Environment/LinearApproximationIntegrator.h
+++ b/Environment/LinearApproximationIntegrator.h
@@ -12,7 +12,7 @@
#include <corsika/geometry/Line.h>
#include <corsika/geometry/Point.h>
-#include <corsika/geometry/Trajectory.h>
+#include <corsika/setup/SetupTrajectory.h>
namespace corsika::environment {
template <class TDerived>
@@ -20,30 +20,28 @@ namespace corsika::environment {
auto const& GetImplementation() const { return *static_cast<TDerived const*>(this); }
public:
- auto IntegrateGrammage(
- corsika::geometry::Trajectory<corsika::geometry::Line> const& line,
- corsika::units::si::LengthType length) const {
+ auto IntegrateGrammage(corsika::setup::Trajectory 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());
+ auto const c1 = GetImplementation().fRho.FirstDerivative(line.GetPosition(0),
+ line.GetDirection(0));
return (c0 + 0.5 * c1 * length) * length;
}
- auto ArclengthFromGrammage(
- corsika::geometry::Trajectory<corsika::geometry::Line> const& line,
- corsika::units::si::GrammageType grammage) const {
+ auto ArclengthFromGrammage(corsika::setup::Trajectory 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());
+ auto const c1 = GetImplementation().fRho.FirstDerivative(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::setup::Trajectory 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/MediumTypes.h b/Environment/MediumTypes.h
new file mode 100644
index 0000000000000000000000000000000000000000..bc44560fd2251927d10f0c2c65d55c60d656d2f6
--- /dev/null
+++ b/Environment/MediumTypes.h
@@ -0,0 +1,22 @@
+/*
+ * (c) Copyright 2020 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
+
+namespace corsika::environment {
+
+ /**
+ * Medium types are useful most importantly for effective models
+ * like energy losses. a particular medium (mixture of components)
+ * may have specif properties not reflected by its mixture of
+ * components.
+ */
+
+ enum class EMediumType { eUnknown, eAir, eWater, eIce, eRock };
+
+} // namespace corsika::environment
diff --git a/Environment/NoMagneticField.h b/Environment/NoMagneticField.h
new file mode 100644
index 0000000000000000000000000000000000000000..74cbd8cfc51edac6db9bd05269f33ddb37aa23b9
--- /dev/null
+++ b/Environment/NoMagneticField.h
@@ -0,0 +1,58 @@
+/*
+ * (c) Copyright 2020 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/IMagneticFieldModel.h>
+#include <corsika/geometry/RootCoordinateSystem.h>
+
+namespace corsika::environment {
+
+ /**
+ * A uniform (constant) magnetic field.
+ *
+ * This class returns the same magnetic field vector
+ * for all evaluated locations.
+ *
+ */
+ template <typename T>
+ class NoMagneticField : public T {
+
+ // a type-alias for a magnetic field vector
+ using MagneticFieldVector =
+ corsika::geometry::Vector<corsika::units::si::magnetic_flux_density_d>;
+
+ public:
+ /**
+ * Construct a NoMagneticField.
+ *
+ * This is initialized with a fixed magnetic field
+ * and returns this magnetic field at all locations.
+ *
+ * @param field The fixed magnetic field to return.
+ */
+ template <typename... Args>
+ NoMagneticField(Args&&... args)
+ : T(std::forward<Args>(args)...) {}
+
+ /**
+ * Evaluate the magnetic field at a given location.
+ *
+ * @param point The location to evaluate the field at.
+ * @returns The magnetic field vector.
+ */
+ MagneticFieldVector GetMagneticField(
+ corsika::geometry::Point const&) const final override {
+ CoordinateSystem const& gCS =
+ RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+ return MagneticFieldVector(gCS, {0_T, 0_T, 0_T});
+ }
+
+ }; // END: class MagneticField
+
+} // namespace corsika::environment
diff --git a/Environment/ShowerAxis.cc b/Environment/ShowerAxis.cc
index 3cc5779e08a6c411a39a7d5aeec0831f5c1e4f2b..76fdc66724ba57850dd404732f0b7940984176fe 100644
--- a/Environment/ShowerAxis.cc
+++ b/Environment/ShowerAxis.cc
@@ -16,14 +16,17 @@ using namespace corsika::units::si;
using namespace corsika;
GrammageType ShowerAxis::X(LengthType l) const {
- auto const fractionalBin = l / steplength_;
+ double const fractionalBin = l / steplength_;
int const lower = fractionalBin; // indices of nearest X support points
- auto const lambda = fractionalBin - lower;
- decltype(X_.size()) const upper = lower + 1;
+ double const frac = fractionalBin - lower;
+ unsigned int const upper = lower + 1;
- if (lower < 0) {
+ if (fractionalBin < 0) {
C8LOG_ERROR("cannot extrapolate to points behind point of injection l={} m", l / 1_m);
- throw std::runtime_error("cannot extrapolate to points behind point of injection");
+ if (throw_) {
+ throw std::runtime_error("cannot extrapolate to points behind point of injection");
+ }
+ return minimumX();
}
if (upper >= X_.size()) {
@@ -31,16 +34,19 @@ GrammageType ShowerAxis::X(LengthType l) const {
fmt::format("shower axis too short, cannot extrapolate (l / max_length_ = {} )",
l / max_length_);
C8LOG_ERROR(err);
- throw std::runtime_error(err.c_str());
+ if (throw_) { throw std::runtime_error(err.c_str()); }
+ return maximumX();
}
- assert(0 <= lambda && lambda <= 1.);
+ C8LOG_TRACE("showerAxis::X frac={}, fractionalBin={}, lower={}, upper={}", frac,
+ fractionalBin, lower, upper);
+ assert(0 <= frac && frac <= 1.);
- C8LOG_TRACE("ShowerAxis::X l={} m, lower={}, lambda={}, upper={}", l / 1_m, lower,
- lambda, upper);
+ C8LOG_TRACE("ShowerAxis::X l={} m, lower={}, frac={}, upper={}", l / 1_m, lower, frac,
+ upper);
// linear interpolation between X[lower] and X[upper]
- return X_[upper] * lambda + X_[lower] * (1 - lambda);
+ return X_[upper] * frac + X_[lower] * (1 - frac);
}
LengthType ShowerAxis::steplength() const { return steplength_; }
diff --git a/Environment/ShowerAxis.h b/Environment/ShowerAxis.h
index bfc6f1212022535a3e727c6ce8afef1a6ee771e6..052aa0bff6ac286039b4fbe580b4093b55d3e396 100644
--- a/Environment/ShowerAxis.h
+++ b/Environment/ShowerAxis.h
@@ -45,9 +45,11 @@ namespace corsika::environment {
template <typename TEnvModel>
ShowerAxis(geometry::Point const& pStart,
geometry::Vector<units::si::length_d> length,
- environment::Environment<TEnvModel> const& env, int steps = 10'000)
+ environment::Environment<TEnvModel> const& env, bool doThrow = false,
+ int steps = 10'000)
: pointStart_(pStart)
, length_(length)
+ , throw_(doThrow)
, max_length_(length_.norm())
, steplength_(max_length_ / steps)
, axis_normalized_(length / max_length_)
@@ -104,6 +106,7 @@ namespace corsika::environment {
private:
geometry::Point const pointStart_;
geometry::Vector<units::si::length_d> const length_;
+ bool throw_ = false;
units::si::LengthType const max_length_, steplength_;
geometry::Vector<units::si::dimensionless_d> const axis_normalized_;
std::vector<units::si::GrammageType> X_;
diff --git a/Environment/SlidingPlanarExponential.h b/Environment/SlidingPlanarExponential.h
index 41058d864f1904b5add45ee902c77ded3c3eb8ec..aae3ccef18ebd1abfafe48afe60963a3ed5da6f1 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,
+ setup::Trajectory 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,
+ setup::Trajectory 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 3acb988e85dbf2f097cf2c3517202af3cf65166f..ad3306810d4efccbee8905d5d281785e6180880d 100644
--- a/Environment/VolumeTreeNode.h
+++ b/Environment/VolumeTreeNode.h
@@ -8,19 +8,18 @@
#pragma once
-#include <corsika/environment/IMediumModel.h>
+#include <corsika/environment/IEmpty.hpp>
#include <corsika/geometry/Volume.h>
#include <memory>
#include <vector>
namespace corsika::environment {
- class Empty {}; //<! intended for usage as default template argument
-
- template <typename TModelProperties = Empty>
+ template <typename TModelProperties = IEmpty>
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>;
@@ -92,7 +91,7 @@ namespace corsika::environment {
fExcludedNodes.push_back(pNode.get());
}
- auto* GetParent() const { return fParentNode; };
+ const VTN_type* GetParent() const { return fParentNode; };
auto const& GetChildNodes() const { return fChildNodes; }
@@ -115,14 +114,15 @@ namespace corsika::environment {
void SetModelProperties(IMPSharedPtr ptr) { fModelProperties = ptr; }
+ /*
template <class MediumType, typename... Args>
static auto CreateMedium(Args&&... args) {
- static_assert(std::is_base_of_v<IMediumModel, MediumType>,
+ static_assert(std::is_base_of_v<, MediumType>,
"unusable type provided, needs to be derived from \"IMediumModel\"");
return std::make_shared<MediumType>(std::forward<Args>(args)...);
}
-
+ */
private:
std::vector<VTNUPtr> fChildNodes;
std::vector<VolumeTreeNode<IModelProperties> const*> fExcludedNodes;
diff --git a/Environment/testEnvironment.cc b/Environment/testEnvironment.cc
index e93eb64f63835c18f50dafc9b2ce8964309e3985..eecae62928a54b4876a164389c5d84626c93cef0 100644
--- a/Environment/testEnvironment.cc
+++ b/Environment/testEnvironment.cc
@@ -28,6 +28,8 @@
#include <corsika/particles/ParticleProperties.h>
#include <corsika/units/PhysicalUnits.h>
+#include <corsika/setup/SetupTrajectory.h>
+
#include <catch2/catch.hpp>
using namespace corsika::geometry;
@@ -61,8 +63,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);
-
+ setup::Trajectory const trajectory = setup::testing::make_track<setup::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 +71,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);
+ setup::Trajectory const trajectory = setup::testing::make_track<setup::Trajectory>(line, tEnd);
LengthType const length = 2 * lambda;
GrammageType const exact = rho0 * lambda * (exp(length / lambda) - 1);
@@ -81,11 +82,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);
+ setup::Trajectory const trajectory = setup::testing::make_track<setup::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 +94,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);
+ setup::Trajectory const trajectory = setup::testing::make_track<setup::Trajectory>(line, tEnd);
double const cosTheta = M_SQRT1_2;
LengthType const length = 2 * lambda;
GrammageType const exact =
@@ -127,7 +128,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);
+ setup::Trajectory const trajectory = setup::testing::make_track<setup::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 +167,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);
+ setup::Trajectory const trajectory = setup::testing::make_track<setup::Trajectory>(line, tEnd);
Exponential const e;
DensityFunction<decltype(e), LinearApproximationIntegrator> const rho(e);
@@ -292,7 +293,7 @@ TEST_CASE("UniformMagneticField w/ Homogeneous Medium") {
auto const tEnd = 1_s;
// and the associated trajectory
- Trajectory<Line> const trajectory(line, tEnd);
+ setup::Trajectory const trajectory = setup::testing::make_track<setup::Trajectory>(line, tEnd);
// and check the integrated grammage
CHECK((medium.IntegratedGrammage(trajectory, 3_m) / (density * 3_m)) == Approx(1));
@@ -395,7 +396,7 @@ TEST_CASE("UniformRefractiveIndex w/ Homogeneous") {
auto const tEnd = 1_s;
// and the associated trajectory
- Trajectory<Line> const trajectory(line, tEnd);
+ setup::Trajectory const trajectory = setup::testing::make_track<setup::Trajectory>(line, tEnd);
// and check the integrated grammage
CHECK((medium.IntegratedGrammage(trajectory, 3_m) / (density * 3_m)) == Approx(1));
@@ -469,7 +470,7 @@ TEST_CASE("MediumPropertyModel w/ Homogeneous") {
auto const tEnd = 1_s;
// and the associated trajectory
- Trajectory<Line> const trajectory(line, tEnd);
+ setup::Trajectory const trajectory = setup::testing::make_track<setup::Trajectory>(line, tEnd);
// and check the integrated grammage
CHECK((medium.IntegratedGrammage(trajectory, 3_m) / (density * 3_m)) == Approx(1));
diff --git a/Environment/testShowerAxis.cc b/Environment/testShowerAxis.cc
index 232b734cb9f4153c73826c258d109d424b83e982..e2862176bc0864bac62238d6b06dae103bccf3b2 100644
--- a/Environment/testShowerAxis.cc
+++ b/Environment/testShowerAxis.cc
@@ -64,7 +64,9 @@ TEST_CASE("Homogeneous Density") {
Point const injectionPos = showerCore + Vector<dimensionless_d>{cs, {0, 0, 1}} * t;
environment::ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos),
- *env, 20};
+ *env,
+ false, // -> do not throw exceptions
+ 20}; // -> number of bins
CHECK(showerAxis.steplength() == 500_m);
diff --git a/Framework/Analytics/testClassTimer.cc b/Framework/Analytics/testClassTimer.cc
index 95f82b4b2caa93394634e5c39ee542dacda94f2f..6185036b5ab36b1a2a203759fa18a06120020bc5 100644
--- a/Framework/Analytics/testClassTimer.cc
+++ b/Framework/Analytics/testClassTimer.cc
@@ -74,7 +74,7 @@ public:
class fooT1 {
public:
template <typename T1, typename T2>
- int inside_t(T1 a, T2 b, T2 c) {
+ int inside_t(T1, T2, T2) {
return 123;
}
};
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index 43a5275435480e45aaa3bb06ac824f2ba5636f40..bc29bb620fea26ae1225f4832395bb6834426fbd 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -19,8 +19,6 @@
#include <corsika/stack/history/EventType.hpp>
#include <corsika/stack/history/HistorySecondaryProducer.hpp>
-#include <corsika/setup/SetupTrajectory.h>
-
/* see Issue 161, we need to include SetupStack only because we need
to globally define StackView. This is clearly not nice and should
be changed, when possible. It might be that StackView needs to be
@@ -33,6 +31,11 @@
#include <cmath>
#include <limits>
+#include <boost/type_index.hpp>
+using boost::typeindex::type_id_with_cvr;
+
+#include <fstream>
+
/**
* The cascade namespace assembles all objects needed to simulate full particles cascades.
*/
@@ -72,17 +75,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;
@@ -104,6 +96,8 @@ namespace corsika::cascade {
}
}
+ ~Cascade(){};
+
/**
* The Run function is the main simulation loop, which processes
* particles from the Stack until the Stack is empty.
@@ -117,7 +111,7 @@ namespace corsika::cascade {
count_++;
auto pNext = stack_.GetNextParticle();
C8LOG_DEBUG(
- "============== next particle : count={}, pid={}, "
+ "============== next particle : count={}, pid={} "
", stack entries={}"
", stack deleted={}",
count_, pNext.GetPID(), stack_.getEntries(), stack_.getDeleted());
@@ -160,10 +154,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);
@@ -182,17 +172,9 @@ namespace corsika::cascade {
// assert that particle stays outside void Universe if it has no
// model properties set
- assert(currentLogicalNode != &*environment_.GetUniverse() ||
- environment_.GetUniverse()->HasModelProperties());
-
- // convert next_step from grammage to length
- LengthType const distance_interact =
- currentLogicalNode->GetModelProperties().ArclengthFromGrammage(step,
- next_interact);
-
- // determine the maximum geometric step length from continuous processes
- LengthType const distance_max = process_sequence_.MaxStepLength(vParticle, step);
- C8LOG_DEBUG("distance_max={} m", distance_max / 1_m);
+ assert((currentLogicalNode != &*environment_.GetUniverse() ||
+ environment_.GetUniverse()->HasModelProperties()) &&
+ "FATAL: The environment model has no valid properties set!");
// determine combined total inverse decay time
InverseTimeType const total_inv_lifetime =
@@ -210,23 +192,43 @@ namespace corsika::cascade {
LengthType const distance_decay = next_decay * vParticle.GetMomentum().norm() /
vParticle.GetEnergy() * units::constants::c;
+ // determine geometric tracking
+ auto [step, nextVol] = tracking_.GetTrack(vParticle);
+ auto geomMaxLength = step.GetLength(1);
+
+ // convert next_step from grammage to length
+ LengthType const distance_interact =
+ currentLogicalNode->GetModelProperties().ArclengthFromGrammage(step,
+ next_interact);
+
+ // determine the maximum geometric step length
+ LengthType const continuous_max_dist = process_sequence_.MaxStepLength(vParticle, step);
+
// take minimum of geometry, interaction, decay for next step
- auto const min_distance =
- std::min({distance_interact, distance_decay, distance_max, geomMaxLength});
+ auto min_distance =
+ std::min({distance_interact, distance_decay, continuous_max_dist, geomMaxLength});
- C8LOG_DEBUG("transport particle by : {} m", min_distance / 1_m);
+ C8LOG_DEBUG(
+ "transport particle by : {} m "
+ "Medium transition after: {} m "
+ "Decay after: {} m "
+ "Interaction after: {} m "
+ "Continuous limit: {} m ",
+ min_distance / 1_m, geomMaxLength / 1_m, distance_decay / 1_m,
+ distance_interact / 1_m, continuous_max_dist / 1_m);
// here the particle is actually moved along the trajectory to new position:
- // std::visit(setup::ParticleUpdate<Particle>{vParticle}, step);
- vParticle.SetPosition(step.PositionFromArclength(min_distance));
- // .... also update time, momentum, direction, ...
- vParticle.SetTime(vParticle.GetTime() + min_distance / units::constants::c);
-
- step.LimitEndTo(min_distance);
+ 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
- if (process_sequence_.DoContinuous(vParticle, step) ==
- process::EProcessReturn::eParticleAbsorbed) {
+ process::EProcessReturn status = process_sequence_.DoContinuous(vParticle, step);
+
+ if (status == process::EProcessReturn::eParticleAbsorbed) {
C8LOG_DEBUG("Cascade: delete absorbed particle PID={} E={} GeV",
vParticle.GetPID(), vParticle.GetEnergy() / 1_GeV);
if (!vParticle.isDeleted()) vParticle.Delete();
@@ -245,7 +247,7 @@ namespace corsika::cascade {
TStackView secondaries(vParticle);
- if (min_distance != distance_max) {
+ if (min_distance < continuous_max_dist) {
/*
Create SecondaryView object on Stack. The data container
remains untouched and identical, and 'projectil' is identical
@@ -258,10 +260,9 @@ namespace corsika::cascade {
[[maybe_unused]] auto projectile = secondaries.GetProjectile();
- if (min_distance == distance_interact) {
+ if (distance_interact < distance_decay) {
interaction(secondaries);
} else {
- assert(min_distance == distance_decay);
decay(secondaries);
// make sure particle actually did decay if it should have done so
if (secondaries.getSize() == 1 &&
@@ -286,20 +287,32 @@ namespace corsika::cascade {
fmt::ptr(numericalNodeAfterStep), fmt::ptr(currentLogicalNode));
return numericalNodeAfterStep == currentLogicalNode;
};
-
- assert(assertion()); // numerical and logical nodes don't match
- } else { // boundary crossing, step is limited by volume boundary
- vParticle.SetNode(nextVol);
- /*
- DoBoundary may delete the particle (or not)
-
- caveat: any changes to vParticle, or even the production
- of new secondaries is currently not passed to ParticleCut,
- thus, particles outside the desired phase space may be produced.
-
- todo: this must be fixed.
- */
- process_sequence_.DoBoundaryCrossing(vParticle, *currentLogicalNode, *nextVol);
+ assert(assertion()); // numerical and logical nodes should
+ // match, we did not cross any volume
+ // boundary
+
+ } else { // boundary crossing, step is limited by volume boundary
+
+ if (nextVol != currentLogicalNode) {
+
+ C8LOG_DEBUG("volume boundary crossing to {}", fmt::ptr(nextVol));
+
+ if (nextVol == environment_.GetUniverse().get()) {
+ C8LOG_DEBUG("particle left physics world, is now in unknown space -> delete");
+ vParticle.Delete();
+ }
+ vParticle.SetNode(nextVol);
+ /*
+ DoBoundary may delete the particle (or not)
+
+ caveat: any changes to vParticle, or even the production
+ of new secondaries is currently not passed to ParticleCut,
+ thus, particles outside the desired phase space may be produced.
+
+ todo: this must be fixed.
+ */
+ process_sequence_.DoBoundaryCrossing(vParticle, *currentLogicalNode, *nextVol);
+ }
}
}
@@ -330,7 +343,7 @@ namespace corsika::cascade {
const auto sample_process = uniDist(rng_);
auto const returnCode = process_sequence_.SelectInteraction(view, sample_process);
if (returnCode != process::EProcessReturn::eInteracted) {
- C8LOG_WARN("Particle did not interace!");
+ C8LOG_WARN("Particle did not interact!");
}
SetEventType(view, history::EventType::Interaction);
return returnCode;
@@ -366,6 +379,17 @@ 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 7cf7b9224327fe3422a86b32e3515f1a44f4af48..a4863d65e68161607b460b87a943382bff677537 100644
--- a/Framework/Cascade/testCascade.cc
+++ b/Framework/Cascade/testCascade.cc
@@ -13,7 +13,6 @@
#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/particles/ParticleProperties.h>
@@ -35,37 +34,66 @@ using namespace corsika::geometry;
#include <limits>
using namespace std;
+/**
+ * testCascade implements an e.m. Heitler model with energy splitting
+ * and a critical energy.
+ *
+ * It resembles one of the most simple cascades you can simulate with CORSIKA8.
+ **/
+
+/*
+ The dummy env (here) doesn't need to have any propoerties
+ */
auto MakeDummyEnv() {
TestEnvironmentType env; // dummy environment
auto& universe = *(env.GetUniverse());
- auto theMedium = TestEnvironmentType::CreateNode<Sphere>(
+ auto world = TestEnvironmentType::CreateNode<Sphere>(
Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
- 100_km * std::numeric_limits<double>::infinity());
+ 1_m * std::numeric_limits<double>::infinity());
- using MyHomogeneousModel = environment::HomogeneousMedium<environment::IMediumModel>;
- theMedium->SetModelProperties<MyHomogeneousModel>(
- 1_g / (1_cm * 1_cm * 1_cm),
- environment::NuclearComposition(
- std::vector<particles::Code>{particles::Code::Proton}, std::vector<float>{1.}));
+ using MyEmptyModel = environment::Empty<environment::IEmpty>;
+ world->SetModelProperties<MyEmptyModel>();
- universe.AddChild(std::move(theMedium));
+ universe.AddChild(std::move(world));
return env;
}
+/**
+ * \class DummyTracking
+ *
+ * For the Heitler model we don't need particle transport.
+ **/
+class DummyTracking {
+
+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;
+ return std::make_tuple(
+ geometry::LineTrajectory(
+ geometry::Line(particle.GetPosition(), initialVelocity),
+ std::numeric_limits<TimeType::value_type>::infinity() * 1_s), // trajectory,
+ // just
+ // go
+ // ahead
+ // forever
+ particle.GetNode()); // next volume node
+ }
+};
+
class ProcessSplit : public process::InteractionProcess<ProcessSplit> {
int fCalls = 0;
- GrammageType fX0;
public:
- ProcessSplit(GrammageType const X0)
- : fX0(X0) {}
-
template <typename Particle>
corsika::units::si::GrammageType GetInteractionLength(Particle const&) const {
- return fX0;
+ return 0_g / square(1_cm);
}
template <typename TSecondaryView>
@@ -73,16 +101,12 @@ public:
fCalls++;
auto const projectile = view.GetProjectile();
const HEPEnergyType E = projectile.GetEnergy();
- view.AddSecondary(
- std::tuple<particles::Code, units::si::HEPEnergyType,
- corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
- projectile.GetPID(), E / 2, projectile.GetMomentum(),
- projectile.GetPosition(), projectile.GetTime()});
- view.AddSecondary(
- std::tuple<particles::Code, units::si::HEPEnergyType,
- corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
- projectile.GetPID(), E / 2, projectile.GetMomentum(),
- projectile.GetPosition(), projectile.GetTime()});
+ view.AddSecondary(std::make_tuple(projectile.GetPID(), E / 2,
+ projectile.GetMomentum(), projectile.GetPosition(),
+ projectile.GetTime()));
+ view.AddSecondary(std::make_tuple(projectile.GetPID(), E / 2,
+ projectile.GetMomentum(), projectile.GetPosition(),
+ projectile.GetTime()));
return EProcessReturn::eInteracted;
}
@@ -122,6 +146,8 @@ public:
TEST_CASE("Cascade", "[Cascade]") {
+ logging::SetLevel(logging::level::trace);
+
HEPEnergyType E0 = 100_GeV;
random::RNGManager& rmng = random::RNGManager::GetInstance();
@@ -129,26 +155,26 @@ TEST_CASE("Cascade", "[Cascade]") {
auto env = MakeDummyEnv();
auto const& rootCS = env.GetCoordinateSystem();
- tracking_line::TrackingLine tracking;
stack_inspector::StackInspector<TestCascadeStack> stackInspect(1, true, E0);
process::NullModel nullModel;
- const GrammageType X0 = 20_g / square(1_cm);
const HEPEnergyType Ecrit = 85_MeV;
- ProcessSplit split(X0);
+ ProcessSplit split;
ProcessCut cut(Ecrit);
auto sequence = process::sequence(nullModel, stackInspect, split, cut);
TestCascadeStack stack;
stack.Clear();
- stack.AddParticle(
- std::tuple<particles::Code, units::si::HEPEnergyType,
- corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
- particles::Code::Electron, E0,
- corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, -1_GeV}),
- Point(rootCS, {0_m, 0_m, 10_km}), 0_ns});
-
- cascade::Cascade<tracking_line::TrackingLine, decltype(sequence), TestCascadeStack,
+ stack.AddParticle(std::make_tuple(
+ particles::Code::Electron, E0,
+ corsika::stack::MomentumVector(
+ rootCS, {0_GeV, 0_GeV,
+ -sqrt(E0 * E0 - units::static_pow<2>(
+ particles::GetMass(particles::Code::Electron)))}),
+ Point(rootCS, {0_m, 0_m, 10_km}), 0_ns));
+
+ DummyTracking tracking;
+ cascade::Cascade<DummyTracking, decltype(sequence), TestCascadeStack,
TestCascadeStackView>
EAS(env, tracking, sequence, stack);
@@ -156,7 +182,7 @@ TEST_CASE("Cascade", "[Cascade]") {
EAS.Run();
CHECK(cut.GetCount() == 2048);
- CHECK(cut.GetCalls() == 2047);
+ CHECK(cut.GetCalls() == 2047); // final particle is still on stack and not yet deleted
CHECK(split.GetCalls() == 2047);
}
diff --git a/Framework/Cascade/testCascade.h b/Framework/Cascade/testCascade.h
index 55740c5601694b2c982c54607bcd4182f519171f..e057441433a58f24e374ee925497fc445ab3cae9 100644
--- a/Framework/Cascade/testCascade.h
+++ b/Framework/Cascade/testCascade.h
@@ -9,14 +9,15 @@
#pragma once
#include <corsika/environment/Environment.h>
+#include <corsika/environment/IEmpty.hpp>
#include <corsika/stack/CombinedStack.h>
#include <corsika/stack/SecondaryView.h>
#include <corsika/stack/node/GeometryNodeStackExtension.h>
#include <corsika/stack/nuclear_extension/NuclearStackExtension.h>
-using TestEnvironmentType =
- corsika::environment::Environment<corsika::environment::IMediumModel>;
+using TestEnvironmentInterface = corsika::environment::IEmpty;
+using TestEnvironmentType = corsika::environment::Environment<TestEnvironmentInterface>;
template <typename T>
using SetupGeometryDataInterface =
diff --git a/Framework/Geometry/CMakeLists.txt b/Framework/Geometry/CMakeLists.txt
index de2f6c5bbd9f69043eaa578a66eb6605faa8e2f9..d302289b0396ed786d3a9c95d3b2aed49145c413 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 01fd0e0f289d978fa8cd69bfaf46031b8a5b2116..5deaa082275d0e459cdd633b421bb3e685742392 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 0000000000000000000000000000000000000000..1d0ca566d8e16a4721619543ac4f024ee50d22ad
--- /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 fd0835308e7f1ca4cb9e0c2192988212eafb0326..f54b02dbd0914b94265b596107037fbac7d0ec43 100644
--- a/Framework/Geometry/Line.h
+++ b/Framework/Geometry/Line.h
@@ -14,6 +14,16 @@
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>;
@@ -22,11 +32,16 @@ namespace corsika::geometry {
VelocityVec const v0;
public:
+ Line() = delete;
+ Line(const Line&) = default;
+ Line(Line&&) = default;
+ Line& operator=(const Line&) = delete;
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 52d3909c39a68bc48105afe829c95a02be9067d0..a2a77efec78e2d6c3b063185c73be8423c265c5a 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,8 @@ 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 1f9c247605de885fc1127cf3bf96843ba8997c85..1580a0db36b6c8e6d973513ac27916bc27a0498e 100644
--- a/Framework/Geometry/Trajectory.h
+++ b/Framework/Geometry/Trajectory.h
@@ -14,42 +14,248 @@
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&) = delete;
+
+ /**
+ * \param theLine The geometric \sa Line object that represents a straight-line
+ * connection
+ *
+ * \param timeLength The time duration to traverse the straight trajectory
+ * in units of \sa TimeType
+ */
+ LineTrajectory(Line const& theLine, corsika::units::si::TimeType timeLength)
+ : line_(theLine)
+ , timeLength_(timeLength)
+ , timeStep_(timeLength)
+ , initialVelocity_(theLine.GetVelocity(corsika::units::si::TimeType::zero()))
+ , finalVelocity_(theLine.GetVelocity(timeLength)) {}
+
+ /**
+ * \param theLine The geometric \sa Line object that represents a straight-line
+ * connection
+ *
+ * \param timeLength The time duration to traverse the straight trajectory
+ * in units of \sa TimeType
+ *
+ * \param timeStep Time duration to folow eventually curved
+ * trajectory in units of \sa TimesType
+ *
+ * \param initialV Initial velocity vector at
+ * start of trajectory \param finalV Final velocity vector at start of trajectory
+ */
+ 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;
+ if (timeStep_ ==
+ std::numeric_limits<corsika::units::si::TimeType::value_type>::infinity() * 1_s)
+ return std::numeric_limits<
+ corsika::units::si::LengthType::value_type>::infinity() *
+ 1_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_s;
+ SetFinalVelocity(GetVelocity(0));
+ timeStep_ = limit;
+ } else {
+ // for infinite steps there can't be a difference between
+ // curved and straight trajectory, this is fundamentally
+ // undefined: assume they are the same (which, i.e. is always correct for a
+ // straight line trajectory).
+ //
+ // Final note: only straight-line trajectories should have
+ // infinite steps! Everything else is ill-defined.
+ if (timeStep_ == std::numeric_limits<
+ corsika::units::si::TimeType::value_type>::infinity() *
+ 1_s ||
+ timeLength_ == std::numeric_limits<
+ corsika::units::si::TimeType::value_type>::infinity() *
+ 1_s) {
+ timeLength_ = limit;
+ timeStep_ = limit;
+ // ...and don't touch velocity
+ } 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_; ///! length of straight step (shortest connecting line)
+ corsika::units::si::TimeType timeStep_; ///! length of bend step (curved)
+ VelocityVec initialVelocity_;
+ VelocityVec finalVelocity_;
+ };
+
+ /**
+ * \class LeapFrogTrajectory
+ *
+ *
+ **/
+
+ class LeapFrogTrajectory {
+
+ using VelocityVec = Vector<corsika::units::si::SpeedType::dimension_type>;
+ typedef corsika::geometry::Vector<corsika::units::si::magnetic_flux_density_d>
+ MagneticFieldVector;
- Trajectory(T const& theT, corsika::units::si::TimeType timeLength)
- : T(theT)
- , fTimeLength(timeLength) {}
+ public:
+ LeapFrogTrajectory() = delete;
+ LeapFrogTrajectory(const LeapFrogTrajectory&) = default;
+ LeapFrogTrajectory(LeapFrogTrajectory&&) = default;
+ LeapFrogTrajectory& operator=(const LeapFrogTrajectory&) = delete;
+ 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(corsika::units::si::TimeType t) const {
- return fTraj.GetPosition(t + fTStart);
- }*/
+ const Line GetLine() const {
+ using namespace corsika::units::si;
+ auto D = GetPosition(1) - GetPosition(0);
+ auto d = D.norm();
+ auto v = initialVelocity_;
+ if (d>1_um) { // if trajectory is ultra-short, we do not
+ // re-calculate velocity, just use initial
+ // value. Otherwise, this is numerically unstable
+ v = D/d * GetVelocity(0).norm();
+ }
+ return Line(GetPosition(0), v);
+ }
+ 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;
+ }
+ VelocityVec GetVelocity(double u) const {
+ return initialVelocity_ +
+ initialVelocity_.cross(magneticfield_) * timeStep_ * u * k_;
+ }
- Point GetPosition(double u) const { return T::GetPosition(fTimeLength * u); }
+ 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 6211b3cc4366a1ea78fcaefbdd326f0078adda0f..b528bf1bf3d4d922092a0e517ba8dd1208184749 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));
}
@@ -262,11 +260,5 @@ TEST_CASE("Trajectories") {
(helix.GetPosition(7_s) - helix.PositionFromArclength(helix.ArcLength(0_s, 7_s)))
.norm()
.magnitude() == Approx(0).margin(absMargin));
-
- auto const t = 1234_s;
- Trajectory<Helix> const base(helix, t);
- CHECK(helix.GetPosition(t).GetCoordinates() == base.GetPosition(1.).GetCoordinates());
-
- CHECK(base.ArcLength(0_s, 1_s) / 1_m == Approx(5));
}
}
diff --git a/Framework/ProcessSequence/ProcessSequence.h b/Framework/ProcessSequence/ProcessSequence.h
index 32d3025940e78be21d15d58f5593f32d4a8f514e..02a3d062ea09792f9c9e4f40d030bd7fc3c02414 100644
--- a/Framework/ProcessSequence/ProcessSequence.h
+++ b/Framework/ProcessSequence/ProcessSequence.h
@@ -215,7 +215,7 @@ namespace corsika::process {
const auto particle = view.parent();
lambda_inv_sum += A_.GetInverseInteractionLength(particle);
// check if we should execute THIS process and then EXIT
- if (lambda_inv_select < lambda_inv_sum) {
+ if (lambda_inv_select <= lambda_inv_sum) {
A_.DoInteraction(view);
return EProcessReturn::eInteracted;
}
@@ -229,7 +229,7 @@ namespace corsika::process {
// if this is not a ContinuousProcess --> evaluate probability
lambda_inv_sum += B_.GetInverseInteractionLength(view.parent());
// check if we should execute THIS process and then EXIT
- if (lambda_inv_select < lambda_inv_sum) {
+ if (lambda_inv_select <= lambda_inv_sum) {
B_.DoInteraction(view);
return EProcessReturn::eInteracted;
}
@@ -279,8 +279,8 @@ namespace corsika::process {
// if this is not a ContinuousProcess --> evaluate probability
decay_inv_sum += A_.GetInverseLifetime(view.parent());
// check if we should execute THIS process and then EXIT
- if (decay_inv_select < decay_inv_sum) { // more pedagogical: rndm_select <
- // decay_inv_sum / decay_inv_tot
+ if (decay_inv_select <= decay_inv_sum) { // more pedagogical: rndm_select <
+ // decay_inv_sum / decay_inv_tot
A_.DoDecay(view);
return EProcessReturn::eDecayed;
}
@@ -294,7 +294,7 @@ namespace corsika::process {
// if this is not a ContinuousProcess --> evaluate probability
decay_inv_sum += B_.GetInverseLifetime(view.parent());
// check if we should execute THIS process and then EXIT
- if (decay_inv_select < decay_inv_sum) {
+ if (decay_inv_select <= decay_inv_sum) {
B_.DoDecay(view);
return EProcessReturn::eDecayed;
}
diff --git a/Framework/Units/PhysicalUnits.h b/Framework/Units/PhysicalUnits.h
index 1f207abbc2f3c08bc7b39c49ca3a5db0719e6e5d..36a835e7e9e072bcce1db01f8980099d864d9afe 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/Framework/Utilities/CMakeLists.txt b/Framework/Utilities/CMakeLists.txt
index 95170d4a35b01b6f44eaf6aff11f47ebf0881e00..3373c01718f97a8676ca0b5c5166710d2861e68f 100644
--- a/Framework/Utilities/CMakeLists.txt
+++ b/Framework/Utilities/CMakeLists.txt
@@ -20,6 +20,7 @@ set (
UTILITIES_SOURCES
COMBoost.cc
CorsikaData.cc
+ quartic.cpp
${CORSIKA_FENV})
set (
@@ -31,6 +32,7 @@ set (
sgn.h
CorsikaFenv.h
MetaProgramming.h
+ quartic.h
)
set (
diff --git a/Framework/Utilities/quartic.cpp b/Framework/Utilities/quartic.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..88da8cc905fc94307d9e89104a9d47add6747702
--- /dev/null
+++ b/Framework/Utilities/quartic.cpp
@@ -0,0 +1,130 @@
+/*
+ * (c) Copyright 2020 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 <cmath>
+#include "quartic.h"
+
+//---------------------------------------------------------------------------
+// solve cubic equation x^3 + a*x^2 + b*x + c
+// x - array of size 3
+// In case 3 real roots: => x[0], x[1], x[2], return 3
+// 2 real roots: x[0], x[1], return 2
+// 1 real root : x[0], x[1] ± i*x[2], return 1
+unsigned int solveP3(double* x, double a, double b, double c) {
+ double a2 = a * a;
+ double q = (a2 - 3 * b) / 9;
+ double r = (a * (2 * a2 - 9 * b) + 27 * c) / 54;
+ double r2 = r * r;
+ double q3 = q * q * q;
+ double A, B;
+ if (r2 < q3) {
+ double t = r / sqrt(q3);
+ if (t < -1) t = -1;
+ if (t > 1) t = 1;
+ t = acos(t);
+ a /= 3;
+ q = -2 * sqrt(q);
+ x[0] = q * cos(t / 3) - a;
+ x[1] = q * cos((t + M_2PI) / 3) - a;
+ x[2] = q * cos((t - M_2PI) / 3) - a;
+ return 3;
+ } else {
+ A = -pow(fabs(r) + sqrt(r2 - q3), 1. / 3);
+ if (r < 0) A = -A;
+ B = (0 == A ? 0 : q / A);
+
+ a /= 3;
+ x[0] = (A + B) - a;
+ x[1] = -0.5 * (A + B) - a;
+ x[2] = 0.5 * sqrt(3.) * (A - B);
+ if (fabs(x[2]) < eps) {
+ x[2] = x[1];
+ return 2;
+ }
+
+ return 1;
+ }
+}
+
+//---------------------------------------------------------------------------
+// solve quartic equation x^4 + a*x^3 + b*x^2 + c*x + d
+// Attention - this function returns dynamically allocated array. It has to be released
+// afterwards.
+DComplex* solve_quartic(double a, double b, double c, double d) {
+ double a3 = -b;
+ double b3 = a * c - 4. * d;
+ double c3 = -a * a * d - c * c + 4. * b * d;
+
+ // cubic resolvent
+ // y^3 − b*y^2 + (ac−4d)*y − a^2*d−c^2+4*b*d = 0
+
+ double x3[3];
+ unsigned int iZeroes = solveP3(x3, a3, b3, c3);
+
+ double q1, q2, p1, p2, D, sqD, y;
+
+ y = x3[0];
+ // The essence - choosing Y with maximal absolute value.
+ if (iZeroes != 1) {
+ if (fabs(x3[1]) > fabs(y)) y = x3[1];
+ if (fabs(x3[2]) > fabs(y)) y = x3[2];
+ }
+
+ // h1+h2 = y && h1*h2 = d <=> h^2 -y*h + d = 0 (h === q)
+
+ D = y * y - 4 * d;
+ if (fabs(D) < eps) // in other words - D==0
+ {
+ q1 = q2 = y * 0.5;
+ // g1+g2 = a && g1+g2 = b-y <=> g^2 - a*g + b-y = 0 (p === g)
+ D = a * a - 4 * (b - y);
+ if (fabs(D) < eps) // in other words - D==0
+ p1 = p2 = a * 0.5;
+
+ else {
+ sqD = sqrt(D);
+ p1 = (a + sqD) * 0.5;
+ p2 = (a - sqD) * 0.5;
+ }
+ } else {
+ sqD = sqrt(D);
+ q1 = (y + sqD) * 0.5;
+ q2 = (y - sqD) * 0.5;
+ // g1+g2 = a && g1*h2 + g2*h1 = c ( && g === p ) Krammer
+ p1 = (a * q1 - c) / (q1 - q2);
+ p2 = (c - a * q2) / (q1 - q2);
+ }
+
+ DComplex* retval = new DComplex[4];
+
+ // solving quadratic eq. - x^2 + p1*x + q1 = 0
+ D = p1 * p1 - 4 * q1;
+ if (D < 0.0) {
+ retval[0].real(-p1 * 0.5);
+ retval[0].imag(sqrt(-D) * 0.5);
+ retval[1] = std::conj(retval[0]);
+ } else {
+ sqD = sqrt(D);
+ retval[0].real((-p1 + sqD) * 0.5);
+ retval[1].real((-p1 - sqD) * 0.5);
+ }
+
+ // solving quadratic eq. - x^2 + p2*x + q2 = 0
+ D = p2 * p2 - 4 * q2;
+ if (D < 0.0) {
+ retval[2].real(-p2 * 0.5);
+ retval[2].imag(sqrt(-D) * 0.5);
+ retval[3] = std::conj(retval[2]);
+ } else {
+ sqD = sqrt(D);
+ retval[2].real((-p2 + sqD) * 0.5);
+ retval[3].real((-p2 - sqD) * 0.5);
+ }
+
+ return retval;
+}
diff --git a/Framework/Utilities/quartic.h b/Framework/Utilities/quartic.h
new file mode 100644
index 0000000000000000000000000000000000000000..ac96afb10e57cff019344d650a1df141ad69b223
--- /dev/null
+++ b/Framework/Utilities/quartic.h
@@ -0,0 +1,66 @@
+/***************************************************************************
+ * Copyright (C) 2016 by Саша Миленковић *
+ * sasa.milenkovic.xyz@gmail.com *
+ * *
+ * This program is free software; you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation; either version 2 of the License, or *
+ * (at your option) any later version. *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * ( http://www.gnu.org/licenses/gpl-3.0.en.html ) *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program; if not, write to the *
+ * Free Software Foundation, Inc., *
+ * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
+ ***************************************************************************/
+
+#ifndef QUARTIC_H_INCLUDED
+#define QUARTIC_H_INCLUDED
+
+#include <complex>
+
+const double PI = 3.141592653589793238463L;
+const double M_2PI = 2 * PI;
+const double eps = 1e-12;
+
+typedef std::complex<double> DComplex;
+
+//---------------------------------------------------------------------------
+// useful for testing
+inline DComplex polinom_2(DComplex x, double a, double b) {
+ // Horner's scheme for x*x + a*x + b
+ return x * (x + a) + b;
+}
+
+//---------------------------------------------------------------------------
+// useful for testing
+inline DComplex polinom_3(DComplex x, double a, double b, double c) {
+ // Horner's scheme for x*x*x + a*x*x + b*x + c;
+ return x * (x * (x + a) + b) + c;
+}
+
+//---------------------------------------------------------------------------
+// useful for testing
+inline DComplex polinom_4(DComplex x, double a, double b, double c, double d) {
+ // Horner's scheme for x*x*x*x + a*x*x*x + b*x*x + c*x + d;
+ return x * (x * (x * (x + a) + b) + c) + d;
+}
+
+//---------------------------------------------------------------------------
+// x - array of size 3
+// In case 3 real roots: => x[0], x[1], x[2], return 3
+// 2 real roots: x[0], x[1], return 2
+// 1 real root : x[0], x[1] ± i*x[2], return 1
+unsigned int solveP3(double* x, double a, double b, double c);
+
+//---------------------------------------------------------------------------
+// solve quartic equation x^4 + a*x^3 + b*x^2 + c*x + d
+// Attention - this function returns dynamically allocated array. It has to be released
+// afterwards.
+DComplex* solve_quartic(double a, double b, double c, double d);
+
+#endif // QUARTIC_H_INCLUDED
diff --git a/Processes/CMakeLists.txt b/Processes/CMakeLists.txt
index 4e9f59df367346d51a9072e6a7bfe95ec87fb4c1..63e66a331dfd12c2f3965933db81979fb1931e8a 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 378cd9259488e8ded32678929e987ba507387e07..4c17d7fb419fa00a55a13b8d21aaa198a11ba05f 100644
--- a/Processes/CONEXSourceCut/testCONEXSourceCut.cc
+++ b/Processes/CONEXSourceCut/testCONEXSourceCut.cc
@@ -67,7 +67,6 @@ TEST_CASE("CONEXSourceCut") {
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;
diff --git a/Processes/NullModel/testNullModel.cc b/Processes/NullModel/testNullModel.cc
new file mode 100644
index 0000000000000000000000000000000000000000..cc142637f1a6517873cf9207107192117e35e2f0
--- /dev/null
+++ b/Processes/NullModel/testNullModel.cc
@@ -0,0 +1,54 @@
+/*
+ * (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 <catch2/catch.hpp>
+
+#include <corsika/process/null_model/NullModel.h>
+
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/RootCoordinateSystem.h>
+#include <corsika/geometry/Vector.h>
+
+#include <corsika/units/PhysicalUnits.h>
+
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+
+using namespace corsika::units::si;
+using namespace corsika::process::null_model;
+using namespace corsika;
+
+TEST_CASE("NullModel", "[processes]") {
+
+ auto [env, csPtr, nodePtr] = setup::testing::setupEnvironment(particles::Code::Oxygen);
+ auto const& cs = *csPtr;
+ [[maybe_unused]] auto const& env_dummy = env;
+ [[maybe_unused]] auto const& node_dummy = nodePtr;
+
+ auto [stack, view] =
+ setup::testing::setupStack(particles::Code::Electron, 0, 0, 100_GeV, nodePtr, cs);
+ [[maybe_unused]] const auto& dummyView = view;
+ auto particle = stack->first();
+
+ geometry::Point const origin(cs, {0_m, 0_m, 0_m});
+ geometry::Vector<units::si::SpeedType::dimension_type> v(cs, 0_m / second, 0_m / second,
+ 1_m / second);
+ geometry::Line line(origin, v);
+ geometry::Trajectory<geometry::Line> track(line, 10_s);
+
+ SECTION("interface") {
+
+ NullModel model(10_m);
+
+ [[maybe_unused]] const process::EProcessReturn ret =
+ model.DoContinuous(particle, track);
+ LengthType const length = model.MaxStepLength(particle, track);
+
+ CHECK((length / 10_m) == Approx(1));
+ }
+}
diff --git a/Processes/ObservationPlane/ObservationPlane.cc b/Processes/ObservationPlane/ObservationPlane.cc
index ab3a8023215528084dd242b656597f74e9989d23..dd8b797436a7bb1c892b75ca742ff6a9e0557ade 100644
--- a/Processes/ObservationPlane/ObservationPlane.cc
+++ b/Processes/ObservationPlane/ObservationPlane.cc
@@ -14,31 +14,40 @@
using namespace corsika::process::observation_plane;
using namespace corsika::units::si;
-ObservationPlane::ObservationPlane(geometry::Plane const& obsPlane,
- std::string const& filename, bool deleteOnHit)
+ObservationPlane::ObservationPlane(
+ geometry::Plane const& obsPlane,
+ geometry::Vector<units::si::dimensionless_d> const& x_axis,
+ std::string const& filename, bool deleteOnHit)
: plane_(obsPlane)
, outputStream_(filename)
, deleteOnHit_(deleteOnHit)
, energy_ground_(0_GeV)
- , count_ground_(0) {
- outputStream_ << "#PDG code, energy / eV, distance to center / m" << std::endl;
+ , count_ground_(0)
+ , xAxis_(x_axis.normalized())
+ , yAxis_(obsPlane.GetNormal().cross(xAxis_)) {
+ outputStream_ << "#PDG code, energy / eV, x distance / m, y distance / m" << std::endl;
}
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;
}
const auto energy = particle.GetEnergy();
+ auto const displacement = trajectory.GetPosition(1) - plane_.GetCenter();
+
outputStream_ << static_cast<int>(particles::GetPDG(particle.GetPID())) << ' '
- << energy / 1_eV << ' '
+ << energy / 1_eV << ' ' << displacement.dot(xAxis_) / 1_m << ' '
+ << displacement.dot(yAxis_) / 1_m
<< (trajectory.GetPosition(1) - plane_.GetCenter()).norm() / 1_m
<< std::endl;
@@ -52,20 +61,78 @@ corsika::process::EProcessReturn ObservationPlane::DoContinuous(
}
}
-LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const&,
+LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const& vParticle,
setup::Trajectory const& trajectory) {
+ int chargeNumber;
+ if (corsika::particles::IsNucleus(vParticle.GetPID())) {
+ chargeNumber = vParticle.GetNuclearZ();
+ } else {
+ chargeNumber = corsika::particles::GetChargeNumber(vParticle.GetPID());
+ }
+ auto const* currentLogicalVolumeNode = vParticle.GetNode();
+ 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.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;
+ } 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;
+ }
+ }
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 dist = (trajectory.GetR0() - pointOfIntersection).norm() * 1.0001;
- C8LOG_TRACE("ObservationPlane::MaxStepLength l={} m", dist / 1_m);
- return dist;
+ auto const pointOfIntersection = trajectory.GetLine().GetPosition(timeOfIntersection);
+ std::cout << " obs plane non b-field " << std::endl;
+ return (trajectory.GetLine().GetR0() - pointOfIntersection).norm() * 1.0001;
}
void ObservationPlane::ShowResults() const {
diff --git a/Processes/ObservationPlane/ObservationPlane.h b/Processes/ObservationPlane/ObservationPlane.h
index 86c5ba0f2100043f5a7c4b9c2ba250b1e3a57a20..3b530f6213a6d88160329c3165e004a3ec6b66b7 100644
--- a/Processes/ObservationPlane/ObservationPlane.h
+++ b/Processes/ObservationPlane/ObservationPlane.h
@@ -26,7 +26,9 @@ namespace corsika::process::observation_plane {
class ObservationPlane : public corsika::process::ContinuousProcess<ObservationPlane> {
public:
- ObservationPlane(geometry::Plane const&, std::string const&, bool = true);
+ ObservationPlane(geometry::Plane const&,
+ geometry::Vector<units::si::dimensionless_d> const&,
+ std::string const&, bool = true);
corsika::process::EProcessReturn DoContinuous(
corsika::setup::Stack::ParticleType& vParticle,
@@ -47,5 +49,6 @@ namespace corsika::process::observation_plane {
units::si::HEPEnergyType energy_ground_ = 0 * units::si::electronvolt;
unsigned int count_ground_ = 0;
+ geometry::Vector<units::si::dimensionless_d> const xAxis_, yAxis_;
};
} // namespace corsika::process::observation_plane
diff --git a/Processes/ObservationPlane/ObservationPlane_interpolation.cc b/Processes/ObservationPlane/ObservationPlane_interpolation.cc
new file mode 100644
index 0000000000000000000000000000000000000000..8460c5f5ba34947adc2e2693de897e69dbbebf7d
--- /dev/null
+++ b/Processes/ObservationPlane/ObservationPlane_interpolation.cc
@@ -0,0 +1,68 @@
+/*
+ * (c) Copyright 2019 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.
+ */
+
+#include <corsika/process/observation_plane/ObservationPlane.h>
+
+#include <fstream>
+
+using namespace corsika::process::observation_plane;
+using namespace corsika::units::si;
+
+ObservationPlane::ObservationPlane(
+ geometry::Plane const& obsPlane,
+ geometry::Vector<units::si::dimensionless_d> const& x_axis,
+ std::string const& filename, bool deleteOnHit)
+ : plane_(obsPlane)
+ , outputStream_(filename)
+ , deleteOnHit_(deleteOnHit)
+ , xAxis_(x_axis.normalized())
+ , yAxis_(obsPlane.GetNormal().cross(xAxis_)) {
+ outputStream_ << "#PDG code, energy / eV, x distance / m, y distance / m" << std::endl;
+}
+
+corsika::process::EProcessReturn ObservationPlane::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());
+
+ if (timeOfIntersection < TimeType::zero()) { return process::EProcessReturn::eOk; }
+
+ if (plane_.IsAbove(trajectory.GetR0()) == plane_.IsAbove(trajectory.GetPosition(1))) {
+ return process::EProcessReturn::eOk;
+ }
+
+ auto const displacement = trajectory.GetPosition(1) - plane_.GetCenter();
+
+ outputStream_ << static_cast<int>(particles::GetPDG(particle.GetPID())) << ' '
+ << particle.GetEnergy() * (1 / 1_eV) << ' '
+ << displacement.dot(xAxis_) / 1_m << ' ' << displacement.dot(yAxis_) / 1_m
+ << ' ' << std::endl;
+
+ if (deleteOnHit_) {
+ return process::EProcessReturn::eParticleAbsorbed;
+ } else {
+ return process::EProcessReturn::eOk;
+ }
+}
+
+LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const&,
+ setup::Trajectory const& trajectory) {
+ TimeType const timeOfIntersection =
+ (plane_.GetCenter() - trajectory.GetR0()).dot(plane_.GetNormal()) /
+ trajectory.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;
+}
diff --git a/Processes/ObservationPlane/testObservationPlane.cc b/Processes/ObservationPlane/testObservationPlane.cc
index 7f3f183a6def750818aa1eb05be4bec81937d587..ddaaaad96e4795df9cb2695fde09a3de24e0c069 100644
--- a/Processes/ObservationPlane/testObservationPlane.cc
+++ b/Processes/ObservationPlane/testObservationPlane.cc
@@ -19,6 +19,9 @@
#include <corsika/particles/ParticleProperties.h>
#include <corsika/units/PhysicalUnits.h>
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+
using namespace corsika::units::si;
using namespace corsika::process::observation_plane;
using namespace corsika;
@@ -27,63 +30,52 @@ using namespace corsika::particles;
TEST_CASE("ContinuousProcess interface", "[proccesses][observation_plane]") {
- auto const& rootCS = RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+ auto [env, csPtr, nodePtr] = setup::testing::setupEnvironment(particles::Code::Oxygen);
+ auto const& cs = *csPtr;
+ [[maybe_unused]] auto const& env_dummy = env;
+ [[maybe_unused]] auto const& node_dummy = nodePtr;
/*
- Test with downward going 1_GeV neutrino, starting at 0,1_m,10m
+ Test with downward going 1_GeV neutrino, starting at 0, 1_m, 10m
ObservationPlane has origin at 0,0,0
*/
- Point const start(rootCS, {0_m, 1_m, 10_m});
- Vector<units::si::SpeedType::dimension_type> vec(rootCS, 0_m / second, 0_m / second,
+ auto [stack, viewPtr] =
+ setup::testing::setupStack(particles::Code::NuE, 0, 0, 1_GeV, nodePtr, cs);
+ [[maybe_unused]] setup::StackView& view = *viewPtr;
+ auto particle = stack->GetNextParticle();
+
+ Point const start(cs, {0_m, 1_m, 10_m});
+ 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);
-
- // setup particle stack, and add primary particle
- setup::Stack stack;
- stack.Clear();
- {
- auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
- return sqrt((Elab - m) * (Elab + m));
- };
- stack.AddParticle(
- std::tuple<Code, units::si::HEPEnergyType, corsika::stack::MomentumVector, Point,
- units::si::TimeType>{
- Code::NuMu, 1_GeV,
- corsika::stack::MomentumVector(
- rootCS, {0_GeV, 0_GeV, -elab2plab(1_GeV, NuMu::GetMass())}),
- Point(rootCS, {1_m, 1_m, 10_m}), 0_ns});
- }
- auto particle = stack.GetNextParticle();
+ setup::Trajectory track =
+ setup::testing::make_track<setup::Trajectory>(line, 12_m / units::constants::c);
+
+ particle.SetPosition(Point(cs, {1_m, 1_m, 10_m})); // moving already along -z
SECTION("horizontal plane") {
- Plane const obsPlane(Point(rootCS, {0_m, 0_m, 0_m}),
- Vector<dimensionless_d>(rootCS, {0., 0., 1.}));
- ObservationPlane obs(obsPlane, "particles.dat", true);
+ Plane const obsPlane(Point(cs, {0_m, 0_m, 0_m}),
+ Vector<dimensionless_d>(cs, {0., 0., 1.}));
+ ObservationPlane obs(obsPlane, Vector<dimensionless_d>(cs, {1., 0., 0.}),
+ "particles.dat", true);
const LengthType length = obs.MaxStepLength(particle, track);
const process::EProcessReturn ret = obs.DoContinuous(particle, track);
REQUIRE(length / 10_m == Approx(1).margin(1e-4));
REQUIRE(ret == process::EProcessReturn::eParticleAbsorbed);
-
- /*
- SECTION("horizontal plane") {
- REQUIRE(true); // todo: we have to check content of output file...
-
- }
- */
}
SECTION("inclined plane") {}
SECTION("transparent plane") {
- Plane const obsPlane(Point(rootCS, {0_m, 0_m, 0_m}),
- Vector<dimensionless_d>(rootCS, {0., 0., 1.}));
- ObservationPlane obs(obsPlane, "particles.dat", false);
+ Plane const obsPlane(Point(cs, {0_m, 0_m, 0_m}),
+ Vector<dimensionless_d>(cs, {0., 0., 1.}));
+ ObservationPlane obs(obsPlane, Vector<dimensionless_d>(cs, {1., 0., 0.}),
+ "particles.dat", false);
const LengthType length = obs.MaxStepLength(particle, track);
const process::EProcessReturn ret = obs.DoContinuous(particle, track);
diff --git a/Processes/ParticleCut/testParticleCut.cc b/Processes/ParticleCut/testParticleCut.cc
index 95e217bc73017b0c39e84b415fa722bd82d287a7..c2343f76373acc2712f553e2803dd4ea10e02644 100644
--- a/Processes/ParticleCut/testParticleCut.cc
+++ b/Processes/ParticleCut/testParticleCut.cc
@@ -170,11 +170,11 @@ TEST_CASE("ParticleCut", "[processes]") {
CHECK(cut.GetCutEnergy() == 0_GeV);
}
- corsika::setup::Trajectory const track{
+ corsika::setup::Trajectory const track = setup::testing::make_track<setup::Trajectory>(
geometry::Line{point0,
geometry::Vector<units::si::SpeedType::dimension_type>{
rootCS, {0_m / second, 0_m / second, -units::constants::c}}},
- 12_m / units::constants::c};
+ 12_m / units::constants::c);
SECTION("cut on DoContinous, just invisibles") {
diff --git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc
index 007be947e5d4c8103da17cabee6b39700120463a..c229e886178806bad80b05764e9e82db1e167cdd 100644
--- a/Processes/Sibyll/testSibyll.cc
+++ b/Processes/Sibyll/testSibyll.cc
@@ -182,7 +182,7 @@ TEST_CASE("SibyllInterface", "[processes]") {
CHECK(pSum.norm() / P0 == Approx(1).margin(0.05));
[[maybe_unused]] const GrammageType length = model.GetInteractionLength(particle);
CHECK(length / 1_g * 1_cm * 1_cm == Approx(88.7).margin(0.1));
- CHECK(view.getSize() == 20);
+ // CHECK(view.getSize() == 20); // also sibyll not stable wrt. to compiler changes
}
SECTION("NuclearInteractionInterface") {
@@ -198,7 +198,7 @@ TEST_CASE("SibyllInterface", "[processes]") {
[[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(view);
[[maybe_unused]] const GrammageType length = model.GetInteractionLength(particle);
CHECK(length / 1_g * 1_cm * 1_cm == Approx(44.2).margin(.1));
- CHECK(view.getSize() == 11);
+ // CHECK(view.getSize() == 11); // also sibyll not stable wrt. to compiler changes
}
SECTION("DecayInterface") {
diff --git a/Processes/StackInspector/testStackInspector.cc b/Processes/StackInspector/testStackInspector.cc
index 9814efa02b7f9848eadb24cca7b0b508ac60a763..3d4cdec2c4b63bfbfb99479848e5711013278930 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/TrackWriter/TrackWriter.cc b/Processes/TrackWriter/TrackWriter.cc
index a3d3edbf39eb27f4aaaf66c6255f6cb879572a9a..2d9825d869f69a1da7956c7631e8ffdb6d85fea2 100644
--- a/Processes/TrackWriter/TrackWriter.cc
+++ b/Processes/TrackWriter/TrackWriter.cc
@@ -28,8 +28,9 @@ namespace corsika::process::track_writer {
using namespace std::string_literals;
fFile.open(fFilename);
- fFile << "# PID, E / eV, start coordinates / m, displacement vector to end / m "s
- << '\n';
+ fFile
+ << "# PID, E / eV, start coordinates / m, displacement vector to end / m, steplength / m "s
+ << '\n';
}
template <>
@@ -47,7 +48,9 @@ namespace corsika::process::track_writer {
<< std::setw(width) << std::scientific << std::setprecision(precision) << start[2] / 1_m
<< std::setw(width) << std::scientific << std::setprecision(precision) << delta[0] / 1_m
<< std::setw(width) << std::scientific << std::setprecision(precision) << delta[1] / 1_m
- << std::setw(width) << std::scientific << std::setprecision(precision) << delta[2] / 1_m << '\n';
+ << std::setw(width) << std::scientific << std::setprecision(precision) << delta[2] / 1_m
+ << std::setw(width) << std::scientific << std::setprecision(precision) << delta.norm() / 1_m
+ << '\n';
// clang-format on
return process::EProcessReturn::eOk;
diff --git a/Processes/Tracking/CMakeLists.txt b/Processes/Tracking/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..088b37a41dbd10cd79e30791a166a92ba5229e39
--- /dev/null
+++ b/Processes/Tracking/CMakeLists.txt
@@ -0,0 +1,44 @@
+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})
+
+# #-- -- -- -- -- -- -- -- -- --
+# #code unit testing
+CORSIKA_ADD_TEST (testTracking)
+target_link_libraries (
+ testTracking
+ ProcessTrackingLine
+ ProcessTrackingLeapFrogStraight
+ ProcessTrackingLeapFrogCurved
+ CORSIKAtesting
+)
diff --git a/Processes/Tracking/Intersect.hpp b/Processes/Tracking/Intersect.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..717ad6ca99647ea1b68e5dce5ad645db0c67449d
--- /dev/null
+++ b/Processes/Tracking/Intersect.hpp
@@ -0,0 +1,137 @@
+/*
+ * (c) Copyright 2020 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/Point.h>
+#include <corsika/geometry/Vector.h>
+#include <corsika/units/PhysicalUnits.h>
+#include <corsika/logging/Logging.h>
+#include <corsika/geometry/Intersections.hpp>
+
+#include <limits>
+
+namespace corsika::process::tracking {
+
+ /**
+ * \class Intersect
+ *
+ * This is a CRTP class to provide a generic volume-tree
+ * intersection for the purpose of tracking.
+ *
+ * It return the closest distance in time to the next geometric
+ * intersection, as well as a pointer to the corresponding new
+ * volume.
+ *
+ * User may provide an optional global step-length limit as
+ * parameter. This may be needd for (simpler) algorithms in magnetic
+ * fields, where tracking errors grow linearly with step-length.
+ * Obviously, in the case of the step-length limit, the returend
+ * "next" volume is just the current one.
+ *
+ **/
+
+ template <typename TDerived>
+ class Intersect {
+
+ protected:
+ template <typename TParticle>
+ auto nextIntersect(
+ const TParticle& particle,
+ corsika::units::si::TimeType step_limit =
+ std::numeric_limits<corsika::units::si::TimeType::value_type>::infinity() *
+ corsika::units::si::second) const {
+ using namespace corsika::units::si;
+ using namespace corsika::geometry;
+
+ typedef
+ typename std::remove_reference<decltype(*particle.GetNode())>::type node_type;
+ node_type& volumeNode =
+ *particle.GetNode(); // current "logical" node, from previous tracking step
+ C8LOG_DEBUG("volumeNode={}, numericallyInside={} ", fmt::ptr(&volumeNode),
+ volumeNode.GetVolume().Contains(particle.GetPosition()));
+
+ // start values:
+ TimeType minTime = step_limit;
+ 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 {}, hasIntersections={} ",
+ fmt::ptr(&volumeNode), fmt::ptr(volumeNode.GetParent()),
+ time_intersections_curr.hasIntersections());
+ if (time_intersections_curr.hasIntersections()) {
+ 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);
+ if (!time_intersections.hasIntersections()) { continue; }
+ 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", AND if
+ // extry point is [[much much]] closer than exit point
+ // (because we might have already numerically "exited" it)!
+ if (t_exit > 0_s && t_entry <= minTime &&
+ -t_entry < t_exit) { // protection agains numerical problem, when we already
+ // _exited_ before
+ // enter chile volume 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);
+ if (!time_intersections.hasIntersections()) { continue; }
+ 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/Tracking/testTracking.cc b/Processes/Tracking/testTracking.cc
new file mode 100644
index 0000000000000000000000000000000000000000..aabd87603fdc888fbc4d5a2a3a39fde50a26adb1
--- /dev/null
+++ b/Processes/Tracking/testTracking.cc
@@ -0,0 +1,155 @@
+/*
+ * (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/process/tracking_leapfrog_straight/Tracking.h>
+#include <corsika/process/tracking_line/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;
+
+typedef corsika::geometry::Vector<corsika::units::si::magnetic_flux_density_d>
+ MagneticFieldVector;
+
+template <typename T>
+int sgn(T val) {
+ return (T(0) < val) - (val < T(0));
+}
+
+/*
+ This is the unified and commond unit test for Tracking:
+
+ - tracking_leapfrog_curved::Tracking
+ - tracking_leapfrog_straight::Tracking
+ - tracking_line::Tracking
+
+ */
+
+TEMPLATE_TEST_CASE("TrackingLeapfrog_Curved", "tracking",
+ tracking_leapfrog_curved::Tracking,
+ tracking_leapfrog_straight::Tracking, tracking_line::Tracking) {
+
+ logging::SetLevel(logging::level::trace);
+
+ const HEPEnergyType P0 = 10_GeV;
+
+ auto PID = GENERATE(as<Code>{}, Code::MuPlus, Code::MuPlus, Code::Gamma);
+ // for algorithms that know magnetic deflections choose: +-50uT, 0uT
+ // otherwise just 0uT
+ auto Bfield = GENERATE(filter(
+ []([[maybe_unused]] MagneticFluxType v) {
+ if constexpr (std::is_same_v<TestType, tracking_line::Tracking>)
+ return v == 0_uT;
+ else
+ return true;
+ },
+ values<MagneticFluxType>({50_uT, 0_uT, -50_uT})));
+ // particle --> (world) --> | --> (target)
+ // true: start inside "world" volume
+ // false: start inside "target" volume
+ 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;
+
+ TestType tracking;
+ 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/CMakeLists.txt b/Processes/TrackingLeapFrogCurved/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..8b3794d695fc1b0ac1e57d98241f7f876ede0833
--- /dev/null
+++ b/Processes/TrackingLeapFrogCurved/CMakeLists.txt
@@ -0,0 +1,37 @@
+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})
+
+# Note: all Tracking Algorithms are tested in testTracking
diff --git a/Processes/TrackingLeapFrogCurved/Tracking.h b/Processes/TrackingLeapFrogCurved/Tracking.h
new file mode 100644
index 0000000000000000000000000000000000000000..da636d6563583096e524a8f08b60150a44169068
--- /dev/null
+++ b/Processes/TrackingLeapFrogCurved/Tracking.h
@@ -0,0 +1,316 @@
+/*
+ * (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
+ * The step is caluculated analytically precisely to reach to the next volume
+ *boundary.
+ **/
+ 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();
+ MagneticFieldVector const& 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 * (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:
+ Tracking()
+ : straightTracking_{tracking_line::Tracking()} {}
+
+ 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();
+
+ // 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
+ MagneticFieldVector const& magneticfield =
+ volumeNode.GetModelProperties().GetMagneticField(position);
+ corsika::units::si::MagneticFluxType const magnitudeB = magneticfield.norm();
+ int const chargeNumber = particle.GetChargeNumber();
+ bool const no_deflection = chargeNumber == 0 || magnitudeB == 0_T;
+
+ if (no_deflection) { return GetLinearTrajectory(particle); }
+
+ HEPMomentumType const pAlongB_delta =
+ (particle.GetMomentum() -
+ particle.GetMomentum().parallelProjectionOnto(magneticfield))
+ .norm();
+
+ if (pAlongB_delta == 0_GeV) {
+ // particle travel along, parallel to magnetic field. Rg is
+ // "0", but for purpose of step limit we return infinity here.
+ C8LOG_TRACE("pAlongB_delta is 0_GeV --> parallel");
+ return GetLinearTrajectory(particle);
+ }
+
+ LengthType const gyroradius =
+ (pAlongB_delta * 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;
+ const TimeType steplimit_time = steplimit / initialVelocity.norm();
+ C8LOG_DEBUG("gyroradius {}, steplimit: {} = {}", gyroradius, steplimit,
+ steplimit_time);
+
+ // traverse the environment volume tree and find next
+ // intersection
+ auto [minTime, minNode] =
+ tracking::Intersect<Tracking>::nextIntersect(particle, steplimit_time);
+
+ 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;
+
+ if (sphere.GetRadius() == 1_km * std::numeric_limits<double>::infinity()) {
+ return geometry::Intersections();
+ }
+
+ const int chargeNumber = particle.GetChargeNumber();
+ const auto& position = particle.GetPosition();
+ MagneticFieldVector const& magneticfield = medium.GetMagneticField(position);
+
+ const geometry::Vector<SpeedType::dimension_type> velocity =
+ particle.GetMomentum() / particle.GetEnergy() * corsika::units::constants::c;
+ const geometry::Vector<dimensionless_d> directionBefore =
+ velocity.normalized(); // determine steplength to next volume
+
+ auto const projectedDirection = directionBefore.cross(magneticfield);
+ auto const projectedDirectionSqrNorm = projectedDirection.GetSquaredNorm();
+ bool const isParallel = (projectedDirectionSqrNorm == 0 * square(1_T));
+
+ if (chargeNumber == 0 || magneticfield.norm() == 0_T || isParallel) {
+ return tracking_line::Tracking::Intersect(particle, sphere, medium);
+ }
+
+ bool const numericallyInside = sphere.Contains(particle.GetPosition());
+
+ const auto absVelocity = velocity.norm();
+ auto energy = particle.GetEnergy();
+ auto k = chargeNumber * corsika::units::constants::cSquared * 1_eV /
+ (absVelocity * energy * 1_V);
+
+ auto const denom =
+ (directionBefore.cross(magneticfield)).GetSquaredNorm() * k * k;
+ const double a =
+ ((directionBefore.cross(magneticfield)).dot(position - sphere.GetCenter()) *
+ k +
+ 1) *
+ 4 / (1_m * 1_m * denom);
+ const double b = directionBefore.dot(position - sphere.GetCenter()) * 8 /
+ (denom * 1_m * 1_m * 1_m);
+ const double c = ((position - sphere.GetCenter()).GetSquaredNorm() -
+ (sphere.GetRadius() * sphere.GetRadius())) *
+ 4 / (denom * 1_m * 1_m * 1_m * 1_m);
+ C8LOG_TRACE("denom={}, a={}, b={}, c={}", denom, a, b, c);
+ std::complex<double>* solutions = solve_quartic(0, a, b, c);
+ LengthType d_enter, d_exit;
+ int first = 0, first_entry = 0, first_exit = 0;
+ for (int i = 0; i < 4; i++) {
+ if (solutions[i].imag() == 0) {
+ LengthType const dist = solutions[i].real() * 1_m;
+ C8LOG_TRACE("Solution (real) for current Volume: {} ", dist);
+ if (numericallyInside) {
+ // there must be an entry (negative) and exit (positive) solution
+ if (dist < -0.0001_m) { // security margin to assure transfer to next
+ // logical volume
+ if (first_entry == 0) {
+ d_enter = dist;
+ } else {
+ d_enter = std::max(d_enter, dist); // closest negative to zero (-1e-4) m
+ }
+ first_entry++;
+
+ } else { // thus, dist >= -0.0001_m
+
+ if (first_exit == 0) {
+ d_exit = dist;
+ } else {
+ d_exit = std::min(d_exit, dist); // closest positive to zero (-1e-4) m
+ }
+ first_exit++;
+ }
+ first = int(first_exit > 0) + int(first_entry > 0);
+
+ } else { // thus, numericallyInside == false
+
+ // both physical solutions (entry, exit) must be positive, and as small as
+ // possible
+ if (dist < -0.0001_m) { // need small numerical margin, to assure transport
+ // into next logical volume
+ continue;
+ }
+ if (first == 0) {
+ d_enter = dist;
+ } else {
+ if (dist < d_enter) {
+ d_exit = d_enter;
+ d_enter = dist;
+ } else {
+ d_exit = dist;
+ }
+ }
+ first++;
+ }
+ } // loop over solutions
+ }
+ delete[] solutions;
+
+ if (first != 2) { // entry and exit points found
+ 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)");
+ }
+
+ protected:
+ /**
+ * Use internally stored class tracking_line::Tracking to
+ * perform a straight line tracking, if no magnetic bendig was
+ * detected.
+ *
+ */
+ template <typename TParticle>
+ auto GetLinearTrajectory(TParticle& particle) {
+
+ using namespace corsika::units::si;
+
+ // perform simple linear tracking
+ auto [straightTrajectory, minNode] = straightTracking_.GetTrack(particle);
+
+ // return as leap-frog trajectory
+ return std::make_tuple(
+ geometry::LeapFrogTrajectory(
+ straightTrajectory.GetLine().GetR0(),
+ straightTrajectory.GetLine().GetV0(),
+ MagneticFieldVector(particle.GetPosition().GetCoordinateSystem(), 0_T,
+ 0_T, 0_T),
+ square(0_m) / (square(1_s) * 1_V),
+ straightTrajectory.GetDuration()), // trajectory
+ minNode); // next volume node
+ }
+
+ protected:
+ tracking_line::Tracking
+ straightTracking_; ///! we want this for neutral and B=0T tracks
+
+ }; // namespace tracking_leapfrog_curved
+
+ } // namespace tracking_leapfrog_curved
+
+} // namespace corsika::process
diff --git a/Processes/TrackingLeapFrogStraight/CMakeLists.txt b/Processes/TrackingLeapFrogStraight/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..fba5345506d54fb40d457f58c5a285fa7119abf4
--- /dev/null
+++ b/Processes/TrackingLeapFrogStraight/CMakeLists.txt
@@ -0,0 +1,36 @@
+set (
+ MODEL_HEADERS
+ Tracking.h
+ )
+
+set (
+ MODEL_NAMESPACE
+ corsika/process/tracking_leapfrog_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})
+
+# Note: all Tracking Algorithms are tested in testTracking
diff --git a/Processes/TrackingLeapFrogStraight/Tracking.h b/Processes/TrackingLeapFrogStraight/Tracking.h
new file mode 100644
index 0000000000000000000000000000000000000000..efa46981fe1894dc8154e77268e086a35bf06943
--- /dev/null
+++ b/Processes/TrackingLeapFrogStraight/Tracking.h
@@ -0,0 +1,235 @@
+/*
+ * (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 <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 straight explicit
+ * `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
+ * (at least) one full leap-frog step to the next volume boundary.
+ *
+ **/
+
+ class Tracking : public tracking_line::Tracking {
+
+ public:
+ /**
+ * \param firstFraction fraction of first leap-frog halve step
+ * relative to full linear step to next volume boundary. This
+ * should not be less than 0.5, otherwise you risk that
+ * particles will never travel from one volume to the next
+ * one. A cross should be possible (even likely). If
+ * firstFraction is too big (~1) the resulting calculated error
+ * will be largest.
+ *
+ */
+ Tracking(double firstFraction = 0.55)
+ : firstFraction_(firstFraction) {}
+
+ 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();
+
+ // check if particle is moving at all
+ const auto absVelocity = initialVelocity.norm();
+ if (absVelocity * 1_s == 0_m) {
+ return std::make_tuple(
+ geometry::LineTrajectory(geometry::Line(initialPosition, initialVelocity),
+ 0_s),
+ volumeNode);
+ }
+
+ // charge of the particle, and magnetic field
+ const int chargeNumber = particle.GetChargeNumber();
+ auto magneticfield =
+ volumeNode->GetModelProperties().GetMagneticField(initialPosition);
+ const auto magnitudeB = magneticfield.GetNorm();
+ C8LOG_DEBUG("field={} uT, chargeNumber={}, magnitudeB={} uT",
+ magneticfield.GetComponents() / 1_uT, chargeNumber, magnitudeB / 1_T);
+ bool const no_deflection = chargeNumber == 0 || magnitudeB == 0_T;
+
+ // 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);
+
+ // if particle is non-deflectable, we are done:
+ if (no_deflection) {
+ C8LOG_DEBUG("no deflection. tracking finished");
+ return std::make_tuple(initialTrack, initialTrackNextVolume);
+ }
+
+ HEPMomentumType const pAlongB_delta =
+ (particle.GetMomentum() -
+ particle.GetMomentum().parallelProjectionOnto(magneticfield))
+ .norm();
+
+ if (pAlongB_delta == 0_GeV) {
+ // particle travel along, parallel to magnetic field. Rg is
+ // "0", but for purpose of step limit we return infinity here.
+ C8LOG_TRACE("pAlongB_delta is 0_GeV --> parallel");
+ return std::make_tuple(initialTrack, initialTrackNextVolume);
+ }
+
+ LengthType const gyroradius =
+ (pAlongB_delta * 1_V /
+ (corsika::units::constants::c * abs(chargeNumber) * magnitudeB * 1_eV));
+
+ // 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
+ 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 geometry::Vector<dimensionless_d> direction = initialVelocity.normalized();
+
+ // avoid any intersections within first halve steplength
+ LengthType const firstHalveSteplength =
+ std::min(steplimit, initialTrackLength * firstFraction_);
+
+ 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)={}, deflection={}",
+ position_mid.GetCoordinates(), new_direction.GetComponents(),
+ new_direction_norm,
+ acos(std::min(1.0, direction.dot(new_direction) / new_direction_norm)) * 180 /
+ M_PI);
+
+ // 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...
+
+ LengthType const finalTrackLength = finalTrack.GetLength(1);
+ LengthType const secondLeapFrogLength = firstHalveSteplength * new_direction_norm;
+
+ // check if volume transition is obvious, OR
+ // for numerical reasons, particles slighly bend "away" from a
+ // volume boundary have a very hard time to cross the border,
+ // thus, if secondLeapFrogLength is just slighly shorter (1e-4m) than
+ // finalTrackLength we better just [extend the
+ // secondLeapFrogLength slightly and] force the volume
+ // crossing:
+ bool const switch_volume = finalTrackLength - 0.0001_m <= secondLeapFrogLength;
+ LengthType const secondHalveStepLength =
+ std::min(secondLeapFrogLength, finalTrackLength);
+
+ C8LOG_DEBUG(
+ "finalTrack(0)={}, finalTrack(1)={}, finalTrackLength={}, "
+ "secondLeapFrogLength={}, secondHalveStepLength={}, "
+ "secondLeapFrogLength-finalTrackLength={}, "
+ "secondHalveStepLength-finalTrackLength={}, "
+ "nextVol={}, transition={}",
+ finalTrack.GetPosition(0).GetCoordinates(),
+ finalTrack.GetPosition(1).GetCoordinates(), finalTrackLength,
+ secondLeapFrogLength, secondHalveStepLength,
+ secondLeapFrogLength - finalTrackLength,
+ secondHalveStepLength - finalTrackLength, fmt::ptr(finalTrackNextVolume),
+ switch_volume);
+
+ // perform the second halve-step
+ auto const new_direction_normalized = new_direction.normalized();
+ const Point finalPosition =
+ position_mid + new_direction_normalized * 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
+ (switch_volume ? finalTrackNextVolume : volumeNode));
+ }
+
+ protected:
+ double firstFraction_;
+ };
+
+ } // namespace tracking_leapfrog_straight
+
+} // namespace corsika::process
diff --git a/Processes/TrackingLine/CMakeLists.txt b/Processes/TrackingLine/CMakeLists.txt
index dd626f8778c5a4139fe042fc651bb5e20eb315fe..75dcaec7e4df4f510adab8f0801413fee77c90e0 100644
--- a/Processes/TrackingLine/CMakeLists.txt
+++ b/Processes/TrackingLine/CMakeLists.txt
@@ -1,11 +1,6 @@
set (
MODEL_HEADERS
- TrackingLine.h
- )
-
-set (
- MODEL_SOURCES
- TrackingLine.cc
+ Tracking.h
)
set (
@@ -13,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
@@ -44,11 +34,4 @@ target_include_directories (
install (FILES ${MODEL_HEADERS} DESTINATION include/${MODEL_NAMESPACE})
-# #-- -- -- -- -- -- -- -- -- --
-# #code unit testing
-CORSIKA_ADD_TEST (testTrackingLine)
-target_link_libraries (
- testTrackingLine
- ProcessTrackingLine
- CORSIKAtesting
-)
+# Note: all Tracking Algorithms are tested in testTracking
diff --git a/Processes/TrackingLine/Tracking.cc b/Processes/TrackingLine/Tracking.cc
new file mode 100644
index 0000000000000000000000000000000000000000..6529bf8b705594809f8ab1515ec4d2891a08aa66
--- /dev/null
+++ b/Processes/TrackingLine/Tracking.cc
@@ -0,0 +1,24 @@
+/*
+ * (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 0000000000000000000000000000000000000000..0dfac370a702e864068fdd523e20ec488886629c
--- /dev/null
+++ b/Processes/TrackingLine/Tracking.h
@@ -0,0 +1,124 @@
+/*
+ * (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/Vector.h>
+#include <corsika/geometry/Intersections.hpp>
+#include <corsika/units/PhysicalUnits.h>
+#include <corsika/logging/Logging.h>
+#include <corsika/process/tracking/Intersect.hpp>
+#include <corsika/geometry/Trajectory.h>
+
+#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&) {
+ 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.cc b/Processes/TrackingLine/TrackingLine.cc
deleted file mode 100644
index 23efd27760e873d6cfd9ba564008bd78a95012af..0000000000000000000000000000000000000000
--- a/Processes/TrackingLine/TrackingLine.cc
+++ /dev/null
@@ -1,60 +0,0 @@
-/*
- * (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/TrackingLine.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 {
-
- std::optional<std::pair<TimeType, TimeType>> TimeOfIntersection(Line const& line,
- Sphere const& sphere) {
- auto const delta = line.GetR0() - sphere.GetCenter();
- auto const v = line.GetV0();
- auto const vSqNorm =
- v.squaredNorm(); // todo: get rid of this by having V0 normalized always
- auto const R = sphere.GetRadius();
-
- auto const vDotDelta = v.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 std::make_pair((-vDotDelta - sqDisc) * invDenom,
- (-vDotDelta + sqDisc) * invDenom);
- } else {
- return {};
- }
- }
-
- TimeType TimeOfIntersection(Line const& vLine, Plane const& vPlane) {
- auto const delta = vPlane.GetCenter() - vLine.GetR0();
- auto const v = vLine.GetV0();
- auto const n = vPlane.GetNormal();
- auto const c = n.dot(v);
-
- if (c.magnitude() == 0) {
- return std::numeric_limits<TimeType::value_type>::infinity() * 1_s;
- } else {
- return n.dot(delta) / c;
- }
- }
-} // namespace corsika::process::tracking_line
diff --git a/Processes/TrackingLine/TrackingLine.h b/Processes/TrackingLine/TrackingLine.h
deleted file mode 100644
index 5e20d920c7e0f875618b06943c244ec0e23d329f..0000000000000000000000000000000000000000
--- a/Processes/TrackingLine/TrackingLine.h
+++ /dev/null
@@ -1,129 +0,0 @@
-/*
- * (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/units/PhysicalUnits.h>
-#include <corsika/logging/Logging.h>
-
-#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() = default;
-
- 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> const velocity =
- p.GetMomentum() / p.GetEnergy() * corsika::units::constants::c;
-
- auto const currentPosition = p.GetPosition();
- C8LOG_DEBUG(
- "TrackingLine pid: {}"
- " , E = {} GeV",
- p.GetPID(), p.GetEnergy() / 1_GeV);
- C8LOG_DEBUG("TrackingLine pos: {}", currentPosition.GetCoordinates());
- C8LOG_DEBUG("TrackingLine E: {} GeV", p.GetEnergy() / 1_GeV);
- C8LOG_DEBUG("TrackingLine p: {} GeV", p.GetMomentum().GetComponents() / 1_GeV);
- C8LOG_DEBUG("TrackingLine v: {} ", velocity.GetComponents());
-
- // to do: include effect of magnetic field
- geometry::Line line(currentPosition, velocity);
-
- auto const* currentLogicalVolumeNode = p.GetNode();
-
- C8LOG_DEBUG("numericallyInside = {} ",
- currentLogicalVolumeNode->GetVolume().Contains(currentPosition));
-
- 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
-
- if (auto opt = TimeOfIntersection(line, sphere); opt.has_value()) {
- auto const [t1, t2] = *opt;
- C8LOG_DEBUG("intersection times: {} s; {} s", t1 / 1_s, t2 / 1_s);
- if (t1.magnitude() > 0)
- intersections.emplace_back(t1, &vtn);
- else if (t2.magnitude() > 0)
- C8LOG_DEBUG("inside other volume");
- }
- };
-
- 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;
- }
-
- C8LOG_DEBUG(" t-intersect: {} ", min);
-
- return std::make_tuple(geometry::Trajectory<geometry::Line>(line, min),
- velocity.norm() * min, minIter->second);
- }
- };
-
- } // namespace tracking_line
-
-} // namespace corsika::process
diff --git a/Processes/TrackingLine/testTrackingLine.cc b/Processes/TrackingLine/testTrackingLine.cc
deleted file mode 100644
index 4cfa2bd603fb52bd6080a1167c80faf82d977b8c..0000000000000000000000000000000000000000
--- a/Processes/TrackingLine/testTrackingLine.cc
+++ /dev/null
@@ -1,98 +0,0 @@
-/*
- * (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_line/TrackingLine.h>
-#include <testTrackingLineStack.h> // test-build, and include file is obtained from CMAKE_CURRENT_SOURCE_DIR
-
-#include <corsika/environment/Environment.h>
-#include <corsika/particles/ParticleProperties.h>
-
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/Sphere.h>
-#include <corsika/geometry/Vector.h>
-
-#include <corsika/setup/SetupTrajectory.h>
-using corsika::setup::Trajectory;
-
-#include <catch2/catch.hpp>
-
-using namespace corsika;
-using namespace corsika::process;
-using namespace corsika::units;
-using namespace corsika::geometry;
-
-#include <iostream>
-using namespace std;
-using namespace corsika::units::si;
-
-TEST_CASE("TrackingLine") {
- environment::Environment<environment::Empty> 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<environment::Empty>::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);
-
- auto const [traj, geomMaxLength, nextVol] = tracking.GetTrack(p);
- [[maybe_unused]] auto dummy_geomMaxLength = geomMaxLength;
- [[maybe_unused]] auto dummy_nextVol = nextVol;
-
- REQUIRE((traj.GetPosition(1.) - Point(cs, 0_m, 0_m, radius))
- .GetComponents(cs)
- .norm()
- .magnitude() == Approx(0).margin(1e-4));
- }
-}
diff --git a/Setup/SetupEnvironment.h b/Setup/SetupEnvironment.h
index f4265e3ddf7342b2972e6317a14bbf946c664ff1..87f0c8e305e70b82d657edb16ed23974368f5855 100644
--- a/Setup/SetupEnvironment.h
+++ b/Setup/SetupEnvironment.h
@@ -13,6 +13,7 @@
#include <corsika/environment/IMediumModel.h>
#include <corsika/environment/IMediumPropertyModel.h>
#include <corsika/environment/IRefractiveIndexModel.h>
+#include <corsika/environment/IMagneticFieldModel.h>
namespace corsika::setup {
@@ -39,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;
@@ -52,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:
@@ -63,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 53fe4d04f7ed20210134d9ae8ecf2a76927a9b9b..2a42db7a649054b5cbe72ae6590174a4f926320a 100644
--- a/Setup/SetupTrajectory.h
+++ b/Setup/SetupTrajectory.h
@@ -12,47 +12,74 @@
#include <corsika/geometry/Line.h>
#include <corsika/geometry/Trajectory.h>
-#include <corsika/units/PhysicalUnits.h>
+#include <corsika/process/tracking_line/Tracking.h>
+#include <corsika/process/tracking_leapfrog_curved/Tracking.h>
+#include <corsika/process/tracking_leapfrog_straight/Tracking.h>
-// #include <variant>
+#include <corsika/units/PhysicalUnits.h>
namespace corsika::setup {
- /// definition of Trajectory base class, to be used in tracking and cascades
- typedef corsika::geometry::Trajectory<corsika::geometry::Line> Trajectory;
+ /**
+ Note/Warning: Tracking and Trajectory must fit together !
- /*
- typedef std::variant<std::monostate, corsika::geometry::Trajectory<Line>,
- corsika::geometry::Trajectory<Helix>>
- Trajectory;
+ tracking_leapfrog_curved::Tracking is the result of the Bachelor
+ thesis of Andre Schmidt, KIT. This is a leap-frog algorithm with
+ an analytical, precise calculation of volume intersections. This
+ algorithm needs a LeapFrogTrajectory.
- /// helper visitor to modify Particle by moving along Trajectory
- template <typename Particle>
- class ParticleUpdate {
+ tracking_leapfrog_straight::Tracking is a more simple and direct
+ leap-frog implementation. The two halve steps are coded explicitly
+ as two straight segments. Intersections with other volumes are
+ calculate only on the straight segments. This algorithm is based
+ on LineTrajectory.
- Particle& fP;
+ tracking_line::Tracking is a pure straight tracker. It is based on
+ LineTrajectory.
+ */
+ typedef corsika::process::tracking_leapfrog_curved::Tracking Tracking;
+ //typedef corsika::process::tracking_leapfrog_straight::Tracking Tracking;
+ //typedef corsika::process::tracking_line::Tracking Tracking;
- public:
- ParticleUpdate(Particle& p)
- : fP(p) {}
- void operator()(std::monostate const&) {}
+ /// definition of Trajectory base class, to be used in tracking and cascades
+ //typedef corsika::geometry::LineTrajectory Trajectory;
+ typedef corsika::geometry::LeapFrogTrajectory Trajectory;
- template <typename T>
- void operator()(T const& trajectory) {
- fP.SetPosition(trajectory.GetPosition(1));
- }
- };
+ /**
+
+ The following section is for unit testing only. Eventually it should
+ be moved to "tests".
+
+
+ */
+
+ namespace testing {
- /// 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 TTrack>
+ TTrack make_track(const corsika::geometry::Line& line,
+ const corsika::units::si::TimeType tEnd);
+
+ template <>
+ inline corsika::geometry::LineTrajectory
+ make_track<corsika::geometry::LineTrajectory>(
+ const corsika::geometry::Line& line, const corsika::units::si::TimeType tEnd) {
+ return corsika::geometry::LineTrajectory(line, tEnd);
}
- template <typename T>
- corsika::units::si::TimeType operator()(T const& trajectory) {
- return trajectory.GetDuration();
+
+ template <>
+ inline corsika::geometry::LeapFrogTrajectory
+ make_track<corsika::geometry::LeapFrogTrajectory>(
+ const corsika::geometry::Line& line, const corsika::units::si::TimeType tEnd) {
+ using namespace corsika::units::si;
+ typedef corsika::geometry::Vector<magnetic_flux_density_d> MagneticFieldVector;
+
+ auto const k = square(0_m) / (square(1_s) * 1_V);
+ return corsika::geometry::LeapFrogTrajectory(
+ line.GetR0(), line.GetV0(),
+ MagneticFieldVector{line.GetR0().GetCoordinateSystem(), 0_T, 0_T, 0_T}, k,
+ tEnd);
}
- };
- */
+
+ } // namespace testing
+
} // namespace corsika::setup
diff --git a/Stack/History/HistoryObservationPlane.cpp b/Stack/History/HistoryObservationPlane.cpp
index b1425f18be3ff20fe68f8ffbaf2f6ff9244d6e78..6a8c94ec7d9e33686383dadf479ef56e413d5522 100644
--- a/Stack/History/HistoryObservationPlane.cpp
+++ b/Stack/History/HistoryObservationPlane.cpp
@@ -28,12 +28,13 @@ 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 +54,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/ThirdParty/CMakeLists.txt b/ThirdParty/CMakeLists.txt
index 19cc00b37ffb1b15d30c7ea97f6b79d554d8b7df..b1553b849d8fbf5d8c5c542b967b82ed663695b3 100644
--- a/ThirdParty/CMakeLists.txt
+++ b/ThirdParty/CMakeLists.txt
@@ -228,14 +228,6 @@ endif (NOT (${USE_CONEX_C8} IN_LIST ThirdPartyChoiceValues))
message (STATUS "USE_CONEX_C8='${USE_CONEX_C8}'")
add_library (C8::ext::conex STATIC IMPORTED GLOBAL)
-# this is from the corsika-data/readLib submodule:
-get_directory_property (Boost_IOSTREAMS_FOUND DIRECTORY ${CORSIKA_DATA}/readLib DEFINITION Boost_IOSTREAMS_FOUND)
-set (conex_boost "")
-set (boost_lib "boost_iostreams")
-if (NOT Boost_IOSTREAMS_FOUND)
- set (conex_boost "NO_BOOST=1")
- set (boost_lib "")
-endif (NOT Boost_IOSTREAMS_FOUND)
if ("x_${USE_CONEX_C8}" STREQUAL "x_SYSTEM")
find_package (CONEX REQUIRED)
@@ -260,7 +252,7 @@ else ()
SOURCE_DIR ${CMAKE_CURRENT_BINARY_DIR}/cxroot/source
INSTALL_DIR ${CMAKE_CURRENT_BINARY_DIR}/cxroot/source
CONFIGURE_COMMAND ""
- BUILD_COMMAND make CX_NO_ROOT=1 CORSIKA_8=1 ${conex_boost} CORSIKA_DATA=${CORSIKA_DATA} all
+ BUILD_COMMAND make CX_NO_ROOT=1 CORSIKA_8=1 CORSIKA_DATA=${CORSIKA_DATA} all
INSTALL_COMMAND ""
BUILD_IN_SOURCE ON
EXCLUDE_FROM_ALL TRUE
@@ -278,26 +270,15 @@ else ()
endif ()
-if (Boost_IOSTREAMS_FOUND)
- set_target_properties (
- C8::ext::conex PROPERTIES
- IMPORTED_LOCATION ${CONEX_PREFIX}/lib/${CMAKE_SYSTEM_NAME}/libCONEXsibyll.a
- IMPORTED_NO_SONAME 1
- SKIP_BUILD_RPATH FALSE
- IMPORTED_LINK_INTERFACE_LIBRARIES boost_iostreams
- INTERFACE_INCLUDE_DIRECTORIES
- $<BUILD_INTERFACE:${CONEX_INCLUDE_DIR}>
- )
-else (Boost_IOSTREAMS_FOUND)
- set_target_properties (
- C8::ext::conex PROPERTIES
- IMPORTED_LOCATION ${CONEX_PREFIX}/lib/${CMAKE_SYSTEM_NAME}/libCONEXsibyll.a
- IMPORTED_NO_SONAME 1
- SKIP_BUILD_RPATH FALSE
- INTERFACE_INCLUDE_DIRECTORIES
- $<BUILD_INTERFACE:${CONEX_INCLUDE_DIR}>
- )
-endif (Boost_IOSTREAMS_FOUND)
+set_target_properties (
+ C8::ext::conex PROPERTIES
+ IMPORTED_LOCATION ${CONEX_PREFIX}/lib/${CMAKE_SYSTEM_NAME}/libCONEXsibyll.a
+ IMPORTED_NO_SONAME 1
+ SKIP_BUILD_RPATH FALSE
+ IMPORTED_LINK_INTERFACE_LIBRARIES bz2
+ INTERFACE_INCLUDE_DIRECTORIES
+ $<BUILD_INTERFACE:${CONEX_INCLUDE_DIR}>
+ )
# libz needed for cnpy, used for SaveHistograms
diff --git a/ThirdParty/cnpy/cnpy.cpp b/ThirdParty/cnpy/cnpy.cpp
index b4c26fccf411b18ae5b11d1c607b436dfdc75225..d2eee53d08820803eb14cc56490d699ede2f68ff 100644
--- a/ThirdParty/cnpy/cnpy.cpp
+++ b/ThirdParty/cnpy/cnpy.cpp
@@ -61,8 +61,8 @@ template<> std::vector<char>& cnpy::operator+=(std::vector<char>& lhs, const cha
void cnpy::parse_npy_header(unsigned char* buffer,size_t& word_size, std::vector<size_t>& shape, bool& fortran_order) {
//std::string magic_string(buffer,6);
- uint8_t major_version = *reinterpret_cast<uint8_t*>(buffer+6);
- uint8_t minor_version = *reinterpret_cast<uint8_t*>(buffer+7);
+ [[maybe_unused]] uint8_t major_version = *reinterpret_cast<uint8_t*>(buffer+6);
+ [[maybe_unused]] uint8_t minor_version = *reinterpret_cast<uint8_t*>(buffer+7);
uint16_t header_len = *reinterpret_cast<uint16_t*>(buffer+8);
std::string header(reinterpret_cast<char*>(buffer+9),header_len);
@@ -196,7 +196,7 @@ cnpy::NpyArray load_the_npz_array(FILE* fp, uint32_t compr_bytes, uint32_t uncom
if(nread != compr_bytes)
throw std::runtime_error("load_the_npy_file: failed fread");
- int err;
+ [[maybe_unused]] int err;
z_stream d_stream;
d_stream.zalloc = Z_NULL;
diff --git a/ThirdParty/phys/units/quantity.hpp b/ThirdParty/phys/units/quantity.hpp
index 970773d598ec4e076467fb990c9fdb610b9c936f..0751923b0ca4eb07a75a0d712523bb69b569b904 100644
--- a/ThirdParty/phys/units/quantity.hpp
+++ b/ThirdParty/phys/units/quantity.hpp
@@ -31,6 +31,7 @@
#include <cmath>
#include <cstdlib>
#include <utility> // std::declval
+#include <type_traits> // std::enable_if
/// namespace phys.
@@ -358,6 +359,10 @@ namespace phys {
static constexpr quantity zero() { return quantity{value_type(0.0)}; }
// static constexpr quantity zero = quantity{ value_type( 0.0 ) };
+ // RU, added conversion to T (often: double) for dimensionless_d
+ template <typename DIM=Dims, std::enable_if_t<std::is_same_v<DIM, dimensionless_d>, int> = 0>
+ operator T() { return m_value; }
+
private:
/**
* private initializing constructor.
diff --git a/Tools/plot_tracks.py b/Tools/plot_tracks.py
new file mode 100644
index 0000000000000000000000000000000000000000..2c9b063b3ec5675019a5934286848a389b6dd79e
--- /dev/null
+++ b/Tools/plot_tracks.py
@@ -0,0 +1,63 @@
+#!/usr/bin/env python3
+
+import os
+import sys, getopt
+import re
+
+# (c) Copyright 2020 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.
+
+# with this script you can plot an animation of output of TrackWriter
+
+tracks_dat="$1"
+if [ -z "$tracks_dat" ]; then
+ echo "usage: $0 <tracks.dat> [output.gif]" >&2
+ exit 1
+fi
+
+directory=$(dirname "$tracks_dat")
+hadrons_dat="$directory/hadrons.dat"
+muons_dat="$directory/muons.dat"
+electrons_dat="$directory/electrons.dat"
+gammas_dat="$directory/gammas.dat"
+
+#if [ "$muons_dat" -ot "$tracks_dat" -o "$hadrons_dat" -ot "$muons_dat" ]; then
+ echo "splitting $tracks_dat into $muons_dat and $hadrons_dat."
+ cat "$tracks_dat" | egrep '^\s+-*13\s' > "$muons_dat"
+ cat "$tracks_dat" | egrep '^\s+-*11\s' > "$electrons_dat"
+ cat "$tracks_dat" | egrep '^\s+-*22\s' > "$gammas_dat"
+ cat "$tracks_dat" | egrep -v '^\s+-*13\s' | egrep -v '^\s+-*11\s' | egrep -v '^\s+-*22\s' > "$hadrons_dat"
+#fi
+
+output="$2"
+if [ -z "$output" ]; then
+ output="$tracks_dat.gif"
+fi
+
+echo "creating $output..."
+
+cat <<EOF | gnuplot
+set term gif animate size 900,900
+set output "$output"
+
+#set zrange [0:40e3]
+#set xrange [-10:10]
+#set yrange [-10:10]
+set xlabel "x / m"
+set ylabel "y / m"
+set zlabel "z / m"
+set title "CORSIKA 8 preliminary"
+
+do for [t=0:359:1] {
+# for separate file per angle:
+# set output sprintf("%03d_$output", t)
+
+ set view 80, t
+ splot "$gammas_dat" u 3:4:5:6:7:8 w vectors nohead lt rgb "orange" t "", "$electrons_dat" u 3:4:5:6:7:8 w vectors nohead lt rgb "blue" t "", "$muons_dat" u 3:4:5:6:7:8 w vectors nohead lt rgb "red" t "", "$hadrons_dat" u 3:4:5:6:7:8 w vectors nohead lc rgb "black" t ""
+}
+EOF
+
+exit $?
diff --git a/do-copyright.py b/do-copyright.py
index 4c2b29cc7a2c79418a5a8517018806e6c0aa1a9e..adc278d1124eafdc943d65937ecdb9f552933a40 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"]
-excludeFiles = ['PhysicalConstants.h','CorsikaFenvOSX.cc', 'sgn.h']
+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