From 2a549e446dc845c5f0a20133c121c53ba6cc22e1 Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Sat, 9 Feb 2019 18:43:12 +0100
Subject: [PATCH] everything, w/o sibyll...
---
Documentation/Examples/CMakeLists.txt | 1 +
Documentation/Examples/cascade_example.cc | 2 +-
Environment/FlatAtmosphere/FlatAtmosphere.h | 2 -
Framework/Cascade/CMakeLists.txt | 8 +-
Framework/Cascade/testCascade.cc | 2 +-
Framework/Geometry/BaseTrajectory.h | 1 -
Framework/Geometry/BaseVector.h | 1 -
Framework/Geometry/CoordinateSystem.cc | 1 -
Framework/Geometry/FourVector.h | 1 -
Framework/Geometry/Helix.h | 1 -
Framework/Geometry/Line.h | 1 -
Framework/Geometry/RootCoordinateSystem.h | 1 -
Framework/Geometry/Sphere.h | 1 -
Framework/Geometry/Trajectory.h | 1 -
Framework/Geometry/Volume.h | 1 -
Framework/Geometry/testFourVector.cc | 1 -
Framework/Particles/ParticleProperties.cc | 1 -
Main/Environment.h | 2 -
Main/Stack.h | 2 -
.../HadronicElasticModel.cc | 189 +++++++++++++++++-
.../HadronicElasticModel.h | 182 ++---------------
Processes/NullModel/NullModel.cc | 20 ++
Processes/NullModel/NullModel.h | 12 +-
Processes/NullModel/testNullModel.cc | 16 +-
Processes/Sibyll/NuclearInteraction.h | 1 -
Processes/Sibyll/nuclib.h | 1 -
Processes/TrackWriter/TrackWriter.cc | 47 ++++-
Processes/TrackWriter/TrackWriter.h | 21 +-
Processes/TrackingLine/CMakeLists.txt | 31 ++-
Processes/TrackingLine/TrackingLine.cc | 136 +++++++++++++
Processes/TrackingLine/TrackingLine.h | 131 ++----------
Processes/TrackingLine/testTrackingLine.cc | 31 +--
.../TrackingLine/testTrackingLineStack.cc | 26 +++
.../TrackingLine/testTrackingLineStack.h | 33 +++
Setup/SetupStack.h | 8 +-
35 files changed, 524 insertions(+), 392 deletions(-)
create mode 100644 Processes/TrackingLine/TrackingLine.cc
create mode 100644 Processes/TrackingLine/testTrackingLineStack.cc
create mode 100644 Processes/TrackingLine/testTrackingLineStack.h
diff --git a/Documentation/Examples/CMakeLists.txt b/Documentation/Examples/CMakeLists.txt
index d1bb8d1d..621cab01 100644
--- a/Documentation/Examples/CMakeLists.txt
+++ b/Documentation/Examples/CMakeLists.txt
@@ -29,6 +29,7 @@ target_link_libraries (cascade_example SuperStupidStack CORSIKAunits CORSIKAlogg
CORSIKAcascade
ProcessStackInspector
ProcessTrackWriter
+ ProcessTrackingLine
ProcessHadronicElasticModel
CORSIKAprocesses
CORSIKAparticles
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 82879ccb..f9d1e95a 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -248,7 +248,7 @@ int main() {
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
// setup processes, decays and interactions
- tracking_line::TrackingLine<setup::Stack> tracking(env);
+ tracking_line::TrackingLine<setup::Stack, setup::Trajectory> tracking(env);
stack_inspector::StackInspector<setup::Stack> p0(true);
random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
diff --git a/Environment/FlatAtmosphere/FlatAtmosphere.h b/Environment/FlatAtmosphere/FlatAtmosphere.h
index 4227e0de..6839def4 100644
--- a/Environment/FlatAtmosphere/FlatAtmosphere.h
+++ b/Environment/FlatAtmosphere/FlatAtmosphere.h
@@ -8,5 +8,3 @@
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
-
-
diff --git a/Framework/Cascade/CMakeLists.txt b/Framework/Cascade/CMakeLists.txt
index fe009b21..64669224 100644
--- a/Framework/Cascade/CMakeLists.txt
+++ b/Framework/Cascade/CMakeLists.txt
@@ -15,13 +15,6 @@ add_library (CORSIKAcascade INTERFACE)
CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAcascade ${CORSIKAcascade_NAMESPACE} ${CORSIKAcascade_HEADERS})
-#target_link_libraries (
-# CORSIKAcascade
-# CORSIKAparticles
-# CORSIKAunits
-# CORSIKAthirdparty # for catch2
-# )
-
# include directive for upstream code
target_include_directories (
CORSIKAcascade
@@ -50,6 +43,7 @@ target_link_libraries (
ProcessSibyll
CORSIKAcascade
ProcessStackInspector
+ ProcessTrackingLine
CORSIKAstackinterface
CORSIKAprocesses
CORSIKAparticles
diff --git a/Framework/Cascade/testCascade.cc b/Framework/Cascade/testCascade.cc
index 91505ad8..1e8a7a4d 100644
--- a/Framework/Cascade/testCascade.cc
+++ b/Framework/Cascade/testCascade.cc
@@ -113,7 +113,7 @@ TEST_CASE("Cascade", "[Cascade]") {
rmng.RegisterRandomStream("cascade");
auto env = MakeDummyEnv();
- tracking_line::TrackingLine<setup::Stack> tracking(env);
+ tracking_line::TrackingLine<setup::Stack, setup::Trajectory> tracking(env);
stack_inspector::StackInspector<setup::Stack> p0(true);
diff --git a/Framework/Geometry/BaseTrajectory.h b/Framework/Geometry/BaseTrajectory.h
index eb04eed3..8d13f97a 100644
--- a/Framework/Geometry/BaseTrajectory.h
+++ b/Framework/Geometry/BaseTrajectory.h
@@ -9,7 +9,6 @@
* the license.
*/
-
#ifndef _include_BASETRAJECTORY_H
#define _include_BASETRAJECTORY_H
diff --git a/Framework/Geometry/BaseVector.h b/Framework/Geometry/BaseVector.h
index 4395fb2a..c41d8edf 100644
--- a/Framework/Geometry/BaseVector.h
+++ b/Framework/Geometry/BaseVector.h
@@ -9,7 +9,6 @@
* the license.
*/
-
#ifndef _include_BASEVECTOR_H_
#define _include_BASEVECTOR_H_
diff --git a/Framework/Geometry/CoordinateSystem.cc b/Framework/Geometry/CoordinateSystem.cc
index 18006260..35e29416 100644
--- a/Framework/Geometry/CoordinateSystem.cc
+++ b/Framework/Geometry/CoordinateSystem.cc
@@ -9,7 +9,6 @@
* the license.
*/
-
#include <corsika/geometry/CoordinateSystem.h>
#include <stdexcept>
diff --git a/Framework/Geometry/FourVector.h b/Framework/Geometry/FourVector.h
index dc4c61e6..0a917a15 100644
--- a/Framework/Geometry/FourVector.h
+++ b/Framework/Geometry/FourVector.h
@@ -9,7 +9,6 @@
* the license.
*/
-
#ifndef _include_corsika_framework_geometry_fourvector_h_
#define _include_corsika_framework_geometry_fourvector_h_
diff --git a/Framework/Geometry/Helix.h b/Framework/Geometry/Helix.h
index 45d1d72a..caf2c34d 100644
--- a/Framework/Geometry/Helix.h
+++ b/Framework/Geometry/Helix.h
@@ -9,7 +9,6 @@
* the license.
*/
-
#ifndef _include_HELIX_H_
#define _include_HELIX_H_
diff --git a/Framework/Geometry/Line.h b/Framework/Geometry/Line.h
index fbd227c1..ce795d8d 100644
--- a/Framework/Geometry/Line.h
+++ b/Framework/Geometry/Line.h
@@ -9,7 +9,6 @@
* the license.
*/
-
#ifndef _include_LINETRAJECTORY_H
#define _include_LINETRAJECTORY_H
diff --git a/Framework/Geometry/RootCoordinateSystem.h b/Framework/Geometry/RootCoordinateSystem.h
index e10690b5..0e641b3a 100644
--- a/Framework/Geometry/RootCoordinateSystem.h
+++ b/Framework/Geometry/RootCoordinateSystem.h
@@ -9,7 +9,6 @@
* the license.
*/
-
#ifndef _include_corsika_geometry_rootcoordinatesystem_h_
#define _include_corsika_geometry_rootcoordinatesystem_h_
diff --git a/Framework/Geometry/Sphere.h b/Framework/Geometry/Sphere.h
index c2fc6b0f..5a38820a 100644
--- a/Framework/Geometry/Sphere.h
+++ b/Framework/Geometry/Sphere.h
@@ -9,7 +9,6 @@
* the license.
*/
-
#ifndef _include_SPHERE_H_
#define _include_SPHERE_H_
diff --git a/Framework/Geometry/Trajectory.h b/Framework/Geometry/Trajectory.h
index 4b2a7f76..5be99e36 100644
--- a/Framework/Geometry/Trajectory.h
+++ b/Framework/Geometry/Trajectory.h
@@ -9,7 +9,6 @@
* the license.
*/
-
#ifndef _include_TRAJECTORY_H
#define _include_TRAJECTORY_H
diff --git a/Framework/Geometry/Volume.h b/Framework/Geometry/Volume.h
index 80f4a572..c6c953d9 100644
--- a/Framework/Geometry/Volume.h
+++ b/Framework/Geometry/Volume.h
@@ -9,7 +9,6 @@
* the license.
*/
-
#ifndef _include_VOLUME_H_
#define _include_VOLUME_H_
diff --git a/Framework/Geometry/testFourVector.cc b/Framework/Geometry/testFourVector.cc
index 9e9b303f..a4063f0e 100644
--- a/Framework/Geometry/testFourVector.cc
+++ b/Framework/Geometry/testFourVector.cc
@@ -9,7 +9,6 @@
* the license.
*/
-
#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
// cpp file
#include <catch2/catch.hpp>
diff --git a/Framework/Particles/ParticleProperties.cc b/Framework/Particles/ParticleProperties.cc
index 266f2e9c..396a10af 100644
--- a/Framework/Particles/ParticleProperties.cc
+++ b/Framework/Particles/ParticleProperties.cc
@@ -9,7 +9,6 @@
* the license.
*/
-
#include <corsika/particles/ParticleProperties.h>
#include <iostream>
diff --git a/Main/Environment.h b/Main/Environment.h
index 4227e0de..6839def4 100644
--- a/Main/Environment.h
+++ b/Main/Environment.h
@@ -8,5 +8,3 @@
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
-
-
diff --git a/Main/Stack.h b/Main/Stack.h
index 4227e0de..6839def4 100644
--- a/Main/Stack.h
+++ b/Main/Stack.h
@@ -8,5 +8,3 @@
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
-
-
diff --git a/Processes/HadronicElasticModel/HadronicElasticModel.cc b/Processes/HadronicElasticModel/HadronicElasticModel.cc
index 943195b6..ab349a6b 100644
--- a/Processes/HadronicElasticModel/HadronicElasticModel.cc
+++ b/Processes/HadronicElasticModel/HadronicElasticModel.cc
@@ -11,15 +11,200 @@
#include <corsika/process/hadronic_elastic_model/HadronicElasticModel.h>
+#include <corsika/environment/Environment.h>
+#include <corsika/environment/NuclearComposition.h>
+#include <corsika/geometry/FourVector.h>
+#include <corsika/random/ExponentialDistribution.h>
+#include <corsika/utl/COMBoost.h>
+
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+
+#include <iomanip>
+#include <iostream>
+
+using namespace corsika::setup;
+using Particle = Stack::ParticleType;
+using Track = Trajectory;
+
namespace corsika::process::HadronicElasticModel {
void HadronicElasticInteraction::Init() {}
HadronicElasticInteraction::HadronicElasticInteraction(
- corsika::environment::Environment const& env,
- corsika::units::si::CrossSectionType x, corsika::units::si::CrossSectionType y)
+ environment::Environment const& env, units::si::CrossSectionType x,
+ units::si::CrossSectionType y)
: fX(x)
, fY(y)
, fEnvironment(env) {}
+ template <>
+ units::si::GrammageType HadronicElasticInteraction::GetInteractionLength(
+ Particle const& p, Track&) {
+ using namespace units::si;
+ if (p.GetPID() == particles::Code::Proton) {
+ auto const* currentNode =
+ fEnvironment.GetUniverse()->GetContainingNode(p.GetPosition());
+ auto const& mediumComposition =
+ currentNode->GetModelProperties().GetNuclearComposition();
+
+ auto const& components = mediumComposition.GetComponents();
+ auto const& fractions = mediumComposition.GetFractions();
+
+ auto const projectileMomentum = p.GetMomentum();
+ auto const projectileMomentumSquaredNorm = projectileMomentum.squaredNorm();
+ auto const projectileEnergy = p.GetEnergy();
+
+ auto const avgCrossSection = [&]() {
+ CrossSectionType avgCrossSection = 0_barn;
+
+ for (size_t i = 0; i < fractions.size(); ++i) {
+ auto const targetMass = particles::GetMass(components[i]);
+ auto const s = detail::static_pow<2>(projectileEnergy + targetMass) -
+ projectileMomentumSquaredNorm;
+ avgCrossSection += CrossSection(s) * fractions[i];
+ }
+
+ std::cout << "avgCrossSection: " << avgCrossSection / 1_mbarn << " mb"
+ << std::endl;
+
+ return avgCrossSection;
+ }();
+
+ auto const avgTargetMassNumber = mediumComposition.GetAverageMassNumber();
+
+ GrammageType const interactionLength =
+ avgTargetMassNumber * units::constants::u / avgCrossSection;
+
+ return interactionLength;
+ } else {
+ return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
+ }
+ }
+
+ template <>
+ process::EProcessReturn HadronicElasticInteraction::DoInteraction(Particle& p, Stack&) {
+ if (p.GetPID() != particles::Code::Proton) { return process::EProcessReturn::eOk; }
+
+ using namespace units::si;
+ using namespace units::constants;
+
+ const auto* currentNode =
+ fEnvironment.GetUniverse()->GetContainingNode(p.GetPosition());
+ const auto& mediumComposition =
+ currentNode->GetModelProperties().GetNuclearComposition();
+ const auto& components = mediumComposition.GetComponents();
+
+ std::vector<units::si::CrossSectionType> cross_section_of_components(
+ mediumComposition.GetComponents().size());
+
+ auto const projectileMomentum = p.GetMomentum();
+ auto const projectileMomentumSquaredNorm = projectileMomentum.squaredNorm();
+ auto const projectileEnergy = p.GetEnergy();
+
+ for (size_t i = 0; i < components.size(); ++i) {
+ auto const targetMass = particles::GetMass(components[i]);
+ auto const s = units::si::detail::static_pow<2>(projectileEnergy + targetMass) -
+ projectileMomentumSquaredNorm;
+ cross_section_of_components[i] = CrossSection(s);
+ }
+
+ const auto targetCode =
+ currentNode->GetModelProperties().SampleTarget(cross_section_of_components, fRNG);
+
+ auto const targetMass = particles::GetMass(targetCode);
+
+ std::uniform_real_distribution phiDist(0., 2 * M_PI);
+
+ geometry::FourVector const projectileLab(projectileEnergy, projectileMomentum);
+ geometry::FourVector const targetLab(
+ targetMass, geometry::Vector<units::si::hepmomentum_d>(
+ projectileMomentum.GetCoordinateSystem(), {0_eV, 0_eV, 0_eV}));
+ utl::COMBoost const boost(projectileLab, targetMass);
+
+ auto const projectileCoM = boost.toCoM(projectileLab);
+ auto const targetCoM = boost.toCoM(targetLab);
+
+ auto const pProjectileCoMSqNorm =
+ projectileCoM.GetSpaceLikeComponents().squaredNorm();
+ auto const pProjectileCoMNorm = sqrt(pProjectileCoMSqNorm);
+
+ auto const eProjectileCoM = projectileCoM.GetTimeLikeComponent();
+ auto const eTargetCoM = targetCoM.GetTimeLikeComponent();
+
+ auto const sqrtS = eProjectileCoM + eTargetCoM;
+ auto const s = units::si::detail::static_pow<2>(sqrtS);
+
+ auto const B = this->B(s);
+ std::cout << B << std::endl;
+
+ random::ExponentialDistribution tDist(1 / B);
+ auto const absT = [&]() {
+ decltype(tDist(fRNG)) absT;
+ auto const maxT = 4 * pProjectileCoMSqNorm;
+
+ do {
+ // |t| cannot become arbitrarily large, max. given by GER eq. (4.16), so we just
+ // throw again until we have an acceptable value. Note that the formula holds in
+ // any frame despite of what is stated in the book.
+ absT = tDist(fRNG);
+ } while (absT >= maxT);
+
+ return absT;
+ }();
+
+ std::cout << "HadronicElasticInteraction: s = " << s * invGeVsq
+ << " GeV²; absT = " << absT * invGeVsq
+ << " GeV² (max./GeV² = " << 4 * invGeVsq * projectileMomentumSquaredNorm
+ << ')' << std::endl;
+
+ auto const theta = 2 * asin(sqrt(absT / (4 * pProjectileCoMSqNorm)));
+ auto const phi = phiDist(fRNG);
+
+ auto const projectileScatteredLab = boost.fromCoM(
+ geometry::FourVector<HEPEnergyType, geometry::Vector<hepmomentum_d>>(
+ eProjectileCoM,
+ geometry::Vector<hepmomentum_d>(projectileMomentum.GetCoordinateSystem(),
+ {pProjectileCoMNorm * sin(theta) * cos(phi),
+ pProjectileCoMNorm * sin(theta) * sin(phi),
+ pProjectileCoMNorm * cos(theta)})));
+
+ p.SetMomentum(projectileScatteredLab.GetSpaceLikeComponents());
+ p.SetEnergy(
+ sqrt(projectileScatteredLab.GetSpaceLikeComponents().squaredNorm() +
+ units::si::detail::static_pow<2>(particles::GetMass(
+ p.GetPID())))); // Don't use energy from boost. It can be smaller than
+ // the momentum due to limited numerical accuracy.
+
+ return process::EProcessReturn::eOk;
+ }
+
+ HadronicElasticInteraction::inveV2 HadronicElasticInteraction::B(eV2 s) const {
+ using namespace units::constants;
+ auto constexpr b_p = 2.3;
+ auto const result =
+ (2 * b_p + 2 * b_p + 4 * pow(s * invGeVsq, gfEpsilon) - 4.2) * invGeVsq;
+ std::cout << "B(" << s << ") = " << result / invGeVsq << " GeV¯²" << std::endl;
+ return result;
+ }
+
+ units::si::CrossSectionType HadronicElasticInteraction::CrossSection(
+ SquaredHEPEnergyType s) const {
+ using namespace units::si;
+ using namespace units::constants;
+ // assuming every target behaves like a proton, fX and fY are universal
+ CrossSectionType const sigmaTotal =
+ fX * pow(s * invGeVsq, gfEpsilon) + fY * pow(s * invGeVsq, -gfEta);
+
+ // according to Schuler & Sjöstrand, PRD 49, 2257 (1994)
+ // (we ignore rho because rho^2 is just ~2 %)
+ auto const sigmaElastic =
+ units::si::detail::static_pow<2>(sigmaTotal) /
+ (16 * M_PI * ConvertHEPToSI<CrossSectionType::dimension_type>(B(s)));
+
+ std::cout << "HEM sigmaTot = " << sigmaTotal / 1_mbarn << " mb" << std::endl;
+ std::cout << "HEM sigmaElastic = " << sigmaElastic / 1_mbarn << " mb" << std::endl;
+ return sigmaElastic;
+ }
+
} // namespace corsika::process::HadronicElasticModel
diff --git a/Processes/HadronicElasticModel/HadronicElasticModel.h b/Processes/HadronicElasticModel/HadronicElasticModel.h
index 8ddc55fa..f0396387 100644
--- a/Processes/HadronicElasticModel/HadronicElasticModel.h
+++ b/Processes/HadronicElasticModel/HadronicElasticModel.h
@@ -12,18 +12,16 @@
#ifndef _include_HadronicElasticInteraction_h
#define _include_HadronicElasticInteraction_h
-#include <corsika/environment/Environment.h>
-#include <corsika/environment/NuclearComposition.h>
-#include <corsika/geometry/FourVector.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/process/InteractionProcess.h>
-#include <corsika/random/ExponentialDistribution.h>
#include <corsika/random/RNGManager.h>
+
#include <corsika/units/PhysicalConstants.h>
#include <corsika/units/PhysicalUnits.h>
-#include <corsika/utl/COMBoost.h>
-#include <iomanip>
+namespace corsika::environment {
+ class Environment;
+}
namespace corsika::process::HadronicElasticModel {
/**
@@ -46,37 +44,16 @@ namespace corsika::process::HadronicElasticModel {
using SquaredHEPEnergyType = decltype(corsika::units::si::HEPEnergyType() *
corsika::units::si::HEPEnergyType());
+ using eV2 = decltype(units::si::detail::static_pow<2>(units::si::electronvolt));
+ using inveV2 = decltype(units::si::detail::static_pow<-2>(units::si::electronvolt));
+
corsika::random::RNG& fRNG =
corsika::random::RNGManager::GetInstance().GetRandomStream(
"HadronicElasticModel");
corsika::environment::Environment const& fEnvironment;
- auto B(decltype(units::si::detail::static_pow<2>(units::si::electronvolt)) s) const {
- using namespace corsika::units::constants;
- auto constexpr b_p = 2.3;
- auto const result =
- (2 * b_p + 2 * b_p + 4 * pow(s * invGeVsq, gfEpsilon) - 4.2) * invGeVsq;
- std::cout << "B(" << s << ") = " << result / invGeVsq << " GeV¯²" << std::endl;
- return result;
- }
-
- corsika::units::si::CrossSectionType CrossSection(SquaredHEPEnergyType s) const {
- using namespace corsika::units::si;
- using namespace corsika::units::constants;
- // assuming every target behaves like a proton, fX and fY are universal
- CrossSectionType const sigmaTotal =
- fX * pow(s * invGeVsq, gfEpsilon) + fY * pow(s * invGeVsq, -gfEta);
-
- // according to Schuler & Sjöstrand, PRD 49, 2257 (1994)
- // (we ignore rho because rho^2 is just ~2 %)
- auto const sigmaElastic =
- units::si::detail::static_pow<2>(sigmaTotal) /
- (16 * M_PI * ConvertHEPToSI<CrossSectionType::dimension_type>(B(s)));
-
- std::cout << "HEM sigmaTot = " << sigmaTotal / 1_mbarn << " mb" << std::endl;
- std::cout << "HEM sigmaElastic = " << sigmaElastic / 1_mbarn << " mb" << std::endl;
- return sigmaElastic;
- }
+ inveV2 B(eV2 s) const;
+ corsika::units::si::CrossSectionType CrossSection(SquaredHEPEnergyType s) const;
public:
HadronicElasticInteraction(corsika::environment::Environment const&,
@@ -86,147 +63,12 @@ namespace corsika::process::HadronicElasticModel {
void Init();
template <typename Particle, typename Track>
- corsika::units::si::GrammageType GetInteractionLength(Particle const& p, Track&) {
- using namespace corsika::units::si;
- if (p.GetPID() == corsika::particles::Code::Proton) {
- auto const* currentNode =
- fEnvironment.GetUniverse()->GetContainingNode(p.GetPosition());
- auto const& mediumComposition =
- currentNode->GetModelProperties().GetNuclearComposition();
-
- auto const& components = mediumComposition.GetComponents();
- auto const& fractions = mediumComposition.GetFractions();
-
- auto const projectileMomentum = p.GetMomentum();
- auto const projectileMomentumSquaredNorm = projectileMomentum.squaredNorm();
- auto const projectileEnergy = p.GetEnergy();
-
- auto const avgCrossSection = [&]() {
- CrossSectionType avgCrossSection = 0_barn;
-
- for (size_t i = 0; i < fractions.size(); ++i) {
- auto const targetMass = corsika::particles::GetMass(components[i]);
- auto const s = detail::static_pow<2>(projectileEnergy + targetMass) -
- projectileMomentumSquaredNorm;
- avgCrossSection += CrossSection(s) * fractions[i];
- }
-
- std::cout << "avgCrossSection: " << avgCrossSection / 1_mbarn << " mb"
- << std::endl;
-
- return avgCrossSection;
- }();
-
- auto const avgTargetMassNumber = mediumComposition.GetAverageMassNumber();
-
- GrammageType const interactionLength =
- avgTargetMassNumber * corsika::units::constants::u / avgCrossSection;
-
- return interactionLength;
- } else {
- return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
- }
- }
+ corsika::units::si::GrammageType GetInteractionLength(Particle const& p, Track&);
template <typename Particle, typename Stack>
- corsika::process::EProcessReturn DoInteraction(Particle& p, Stack&) {
- if (p.GetPID() != corsika::particles::Code::Proton) {
- return process::EProcessReturn::eOk;
- }
-
- using namespace units::si;
- using namespace units::constants;
-
- const auto* currentNode =
- fEnvironment.GetUniverse()->GetContainingNode(p.GetPosition());
- const auto& mediumComposition =
- currentNode->GetModelProperties().GetNuclearComposition();
- const auto& components = mediumComposition.GetComponents();
-
- std::vector<units::si::CrossSectionType> cross_section_of_components(
- mediumComposition.GetComponents().size());
-
- auto const projectileMomentum = p.GetMomentum();
- auto const projectileMomentumSquaredNorm = projectileMomentum.squaredNorm();
- auto const projectileEnergy = p.GetEnergy();
-
- for (size_t i = 0; i < components.size(); ++i) {
- auto const targetMass = corsika::particles::GetMass(components[i]);
- auto const s = units::si::detail::static_pow<2>(projectileEnergy + targetMass) -
- projectileMomentumSquaredNorm;
- cross_section_of_components[i] = CrossSection(s);
- }
-
- const auto targetCode = currentNode->GetModelProperties().SampleTarget(
- cross_section_of_components, fRNG);
-
- auto const targetMass = corsika::particles::GetMass(targetCode);
-
- std::uniform_real_distribution phiDist(0., 2 * M_PI);
-
- geometry::FourVector const projectileLab(projectileEnergy, projectileMomentum);
- geometry::FourVector const targetLab(
- targetMass, geometry::Vector<units::si::hepmomentum_d>(
- projectileMomentum.GetCoordinateSystem(), {0_eV, 0_eV, 0_eV}));
- utl::COMBoost const boost(projectileLab, targetMass);
-
- auto const projectileCoM = boost.toCoM(projectileLab);
- auto const targetCoM = boost.toCoM(targetLab);
-
- auto const pProjectileCoMSqNorm =
- projectileCoM.GetSpaceLikeComponents().squaredNorm();
- auto const pProjectileCoMNorm = sqrt(pProjectileCoMSqNorm);
-
- auto const eProjectileCoM = projectileCoM.GetTimeLikeComponent();
- auto const eTargetCoM = targetCoM.GetTimeLikeComponent();
-
- auto const sqrtS = eProjectileCoM + eTargetCoM;
- auto const s = units::si::detail::static_pow<2>(sqrtS);
-
- auto const B = this->B(s);
- std::cout << B << std::endl;
-
- corsika::random::ExponentialDistribution tDist(1 / B);
- auto const absT = [&]() {
- decltype(tDist(fRNG)) absT;
- auto const maxT = 4 * pProjectileCoMSqNorm;
-
- do {
- // |t| cannot become arbitrarily large, max. given by GER eq. (4.16), so we just
- // throw again until we have an acceptable value. Note that the formula holds in
- // any frame despite of what is stated in the book.
- absT = tDist(fRNG);
- } while (absT >= maxT);
-
- return absT;
- }();
-
- std::cout << "HadronicElasticInteraction: s = " << s * invGeVsq
- << " GeV²; absT = " << absT * invGeVsq
- << " GeV² (max./GeV² = " << 4 * invGeVsq * projectileMomentumSquaredNorm
- << ')' << std::endl;
-
- auto const theta = 2 * asin(sqrt(absT / (4 * pProjectileCoMSqNorm)));
- auto const phi = phiDist(fRNG);
-
- auto const projectileScatteredLab = boost.fromCoM(
- geometry::FourVector<HEPEnergyType, geometry::Vector<hepmomentum_d>>(
- eProjectileCoM,
- geometry::Vector<hepmomentum_d>(projectileMomentum.GetCoordinateSystem(),
- {pProjectileCoMNorm * sin(theta) * cos(phi),
- pProjectileCoMNorm * sin(theta) * sin(phi),
- pProjectileCoMNorm * cos(theta)})));
-
- p.SetMomentum(projectileScatteredLab.GetSpaceLikeComponents());
- p.SetEnergy(
- sqrt(projectileScatteredLab.GetSpaceLikeComponents().squaredNorm() +
- units::si::detail::static_pow<2>(corsika::particles::GetMass(
- p.GetPID())))); // Don't use energy from boost. It can be smaller than
- // the momentum due to limited numerical accuracy.
-
- return process::EProcessReturn::eOk;
- }
+ corsika::process::EProcessReturn DoInteraction(Particle& p, Stack&);
};
+
} // namespace corsika::process::HadronicElasticModel
#endif
diff --git a/Processes/NullModel/NullModel.cc b/Processes/NullModel/NullModel.cc
index 09956e7c..c99a464b 100644
--- a/Processes/NullModel/NullModel.cc
+++ b/Processes/NullModel/NullModel.cc
@@ -10,5 +10,25 @@
*/
#include <corsika/process/null_model/NullModel.h>
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+
+using namespace corsika;
+using namespace corsika::process::null_model;
void corsika::process::null_model::NullModel::Init() {}
+
+NullModel::NullModel(units::si::LengthType maxStepLength)
+ : fMaxStepLength(maxStepLength) {}
+
+template <>
+process::EProcessReturn NullModel::DoContinuous(setup::Stack::ParticleType&,
+ setup::Trajectory&, setup::Stack&) const {
+ return process::EProcessReturn::eOk;
+}
+
+template <>
+units::si::LengthType NullModel::MaxStepLength(setup::Stack::ParticleType&,
+ setup::Trajectory&) const {
+ return fMaxStepLength;
+}
diff --git a/Processes/NullModel/NullModel.h b/Processes/NullModel/NullModel.h
index 6222e290..5579d69f 100644
--- a/Processes/NullModel/NullModel.h
+++ b/Processes/NullModel/NullModel.h
@@ -13,7 +13,6 @@
#define _Physics_NullModel_NullModel_h_
#include <corsika/process/ContinuousProcess.h>
-#include <corsika/setup/SetupTrajectory.h>
namespace corsika::process::null_model {
@@ -22,20 +21,15 @@ namespace corsika::process::null_model {
public:
NullModel(corsika::units::si::LengthType maxStepLength =
- corsika::units::si::meter * std::numeric_limits<double>::infinity())
- : fMaxStepLength(maxStepLength) {}
+ corsika::units::si::meter * std::numeric_limits<double>::infinity());
void Init();
template <typename Particle, typename Track, typename Stack>
- process::EProcessReturn DoContinuous(Particle&, Track&, Stack&) const {
- return EProcessReturn::eOk;
- }
+ process::EProcessReturn DoContinuous(Particle&, Track&, Stack&) const;
template <typename Particle, typename Track>
- corsika::units::si::LengthType MaxStepLength(Particle&, Track&) const {
- return fMaxStepLength;
- }
+ corsika::units::si::LengthType MaxStepLength(Particle&, Track&) const;
};
} // namespace corsika::process::null_model
diff --git a/Processes/NullModel/testNullModel.cc b/Processes/NullModel/testNullModel.cc
index 9be839d9..ed639eca 100644
--- a/Processes/NullModel/testNullModel.cc
+++ b/Processes/NullModel/testNullModel.cc
@@ -39,13 +39,15 @@ TEST_CASE("NullModel", "[processes]") {
geometry::Trajectory<geometry::Line> track(line, 10_s);
setup::Stack stack;
- auto particle = stack.AddParticle(
- std::tuple<corsika::particles::Code, corsika::units::si::HEPEnergyType,
- corsika::stack::MomentumVector, corsika::geometry::Point,
- corsika::units::si::TimeType>{
- particles::Code::Electron, 1.5_GeV,
- stack::MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
- geometry::Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s});
+ setup::Stack::ParticleType
+ // auto
+ particle = stack.AddParticle(
+ std::tuple<corsika::particles::Code, corsika::units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, corsika::geometry::Point,
+ corsika::units::si::TimeType>{
+ particles::Code::Electron, 1.5_GeV,
+ stack::MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+ geometry::Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s});
SECTION("interface") {
diff --git a/Processes/Sibyll/NuclearInteraction.h b/Processes/Sibyll/NuclearInteraction.h
index d9eca396..cd273502 100644
--- a/Processes/Sibyll/NuclearInteraction.h
+++ b/Processes/Sibyll/NuclearInteraction.h
@@ -9,7 +9,6 @@
* the license.
*/
-
#ifndef _corsika_process_sibyll_nuclearinteraction_h_
#define _corsika_process_sibyll_nuclearinteraction_h_
diff --git a/Processes/Sibyll/nuclib.h b/Processes/Sibyll/nuclib.h
index 86a2a8f8..e4b89b52 100644
--- a/Processes/Sibyll/nuclib.h
+++ b/Processes/Sibyll/nuclib.h
@@ -9,7 +9,6 @@
* the license.
*/
-
#ifndef _include_nuclib_interface_h_
#define _include_nuclib_interface_h_
diff --git a/Processes/TrackWriter/TrackWriter.cc b/Processes/TrackWriter/TrackWriter.cc
index 65dec50c..b58a40fb 100644
--- a/Processes/TrackWriter/TrackWriter.cc
+++ b/Processes/TrackWriter/TrackWriter.cc
@@ -10,12 +10,45 @@
*/
#include <corsika/process/track_writer/TrackWriter.h>
-#include <string>
-void corsika::process::TrackWriter::TrackWriter::Init() {
- using namespace std::string_literals;
+#include <corsika/particles/ParticleProperties.h>
- fFile.open(fFilename);
- fFile << "# PID, E / eV, start coordinates / m, displacement vector to end / m "s
- << '\n';
-}
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+
+#include <limits>
+
+using namespace corsika::setup;
+using Particle = Stack::ParticleType;
+using Track = Trajectory;
+
+namespace corsika::process::TrackWriter {
+
+ void TrackWriter::Init() {
+ using namespace std::string_literals;
+
+ fFile.open(fFilename);
+ fFile << "# PID, E / eV, start coordinates / m, displacement vector to end / m "s
+ << '\n';
+ }
+
+ template <>
+ process::EProcessReturn TrackWriter::DoContinuous(Particle& p, Track& t, Stack&) {
+ using namespace units::si;
+ auto const start = t.GetPosition(0).GetCoordinates();
+ auto const delta = t.GetPosition(1).GetCoordinates() - start;
+ auto const& name = particles::GetName(p.GetPID());
+
+ fFile << name << " " << p.GetEnergy() / 1_eV << ' ' << start[0] / 1_m << ' '
+ << start[1] / 1_m << ' ' << start[2] / 1_m << " " << delta[0] / 1_m << ' '
+ << delta[1] / 1_m << ' ' << delta[2] / 1_m << '\n';
+
+ return process::EProcessReturn::eOk;
+ }
+
+ template <>
+ units::si::LengthType TrackWriter::MaxStepLength(Particle&, Track&) {
+ return units::si::meter * std::numeric_limits<double>::infinity();
+ }
+
+} // namespace corsika::process::TrackWriter
diff --git a/Processes/TrackWriter/TrackWriter.h b/Processes/TrackWriter/TrackWriter.h
index c61d2a6c..fa243056 100644
--- a/Processes/TrackWriter/TrackWriter.h
+++ b/Processes/TrackWriter/TrackWriter.h
@@ -12,12 +12,10 @@
#ifndef _Processes_TrackWriter_h_
#define _Processes_TrackWriter_h_
-#include <corsika/particles/ParticleProperties.h>
#include <corsika/process/ContinuousProcess.h>
-#include <corsika/setup/SetupTrajectory.h>
#include <corsika/units/PhysicalUnits.h>
+
#include <fstream>
-#include <limits>
#include <string>
namespace corsika::process::TrackWriter {
@@ -31,23 +29,10 @@ namespace corsika::process::TrackWriter {
void Init();
template <typename Particle, typename Track, typename Stack>
- corsika::process::EProcessReturn DoContinuous(Particle& p, Track& t, Stack&) {
- using namespace corsika::units::si;
- auto const start = t.GetPosition(0).GetCoordinates();
- auto const delta = t.GetPosition(1).GetCoordinates() - start;
- auto const& name = corsika::particles::GetName(p.GetPID());
-
- fFile << name << " " << p.GetEnergy() / 1_eV << ' ' << start[0] / 1_m << ' '
- << start[1] / 1_m << ' ' << start[2] / 1_m << " " << delta[0] / 1_m << ' '
- << delta[1] / 1_m << ' ' << delta[2] / 1_m << '\n';
-
- return corsika::process::EProcessReturn::eOk;
- }
+ corsika::process::EProcessReturn DoContinuous(Particle& p, Track& t, Stack&);
template <typename Particle, typename Track>
- corsika::units::si::LengthType MaxStepLength(Particle&, Track&) {
- return corsika::units::si::meter * std::numeric_limits<double>::infinity();
- }
+ corsika::units::si::LengthType MaxStepLength(Particle&, Track&);
private:
std::string const fFilename;
diff --git a/Processes/TrackingLine/CMakeLists.txt b/Processes/TrackingLine/CMakeLists.txt
index ce19ed67..c8b95053 100644
--- a/Processes/TrackingLine/CMakeLists.txt
+++ b/Processes/TrackingLine/CMakeLists.txt
@@ -4,14 +4,37 @@ set (
TrackingLine.h
)
+set (
+ MODEL_SOURCES
+ TrackingLine.cc
+ )
+
set (
MODEL_NAMESPACE
corsika/process/tracking_line
)
-add_library (ProcessTrackingLine INTERFACE)
+add_library (ProcessTrackingLine STATIC ${MODEL_SOURCES})
CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessTrackingLine ${MODEL_NAMESPACE} ${MODEL_HEADERS})
+set_target_properties (
+ ProcessTrackingLine
+ PROPERTIES
+ VERSION ${PROJECT_VERSION}
+ SOVERSION 1
+# PUBLIC_HEADER "${MODEL_HEADERS}"
+ )
+
+# target dependencies on other libraries (also the header onlys)
+target_link_libraries (
+ ProcessTrackingLine
+ CORSIKAsetup
+ CORSIKAutilities
+ CORSIKAenvironment
+ CORSIKAunits
+ CORSIKAenvironment
+ CORSIKAgeometry
+ )
target_include_directories (
ProcessTrackingLine
@@ -29,12 +52,6 @@ add_executable (testTrackingLine testTrackingLine.cc)
target_link_libraries (
testTrackingLine
ProcessTrackingLine
- CORSIKAsetup
- CORSIKAutilities
- CORSIKAenvironment
- CORSIKAunits
- CORSIKAenvironment
- CORSIKAgeometry
CORSIKAthirdparty # for catch2
)
CORSIKA_ADD_TEST(testTrackingLine)
diff --git a/Processes/TrackingLine/TrackingLine.cc b/Processes/TrackingLine/TrackingLine.cc
new file mode 100644
index 00000000..4af036cc
--- /dev/null
+++ b/Processes/TrackingLine/TrackingLine.cc
@@ -0,0 +1,136 @@
+#include <corsika/process/tracking_line/TrackingLine.h>
+
+#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 <algorithm>
+#include <iostream>
+#include <stdexcept>
+#include <utility>
+
+using namespace corsika;
+
+namespace corsika::process::tracking_line {
+
+ template <class Stack, class Trajectory>
+ std::optional<std::pair<corsika::units::si::TimeType, corsika::units::si::TimeType>>
+ TrackingLine<Stack, Trajectory>::TimeOfIntersection(corsika::geometry::Line const& line,
+ geometry::Sphere const& sphere) {
+ using namespace corsika::units::si;
+ auto const& cs = fEnvironment.GetCoordinateSystem();
+ geometry::Point const origin(cs, 0_m, 0_m, 0_m);
+
+ auto const r0 = (line.GetR0() - origin);
+ auto const v0 = line.GetV0();
+ auto const c0 = (sphere.GetCenter() - origin);
+
+ auto const alpha = r0.dot(v0) - 2 * v0.dot(c0);
+ auto const beta = c0.squaredNorm() + r0.squaredNorm() + 2 * c0.dot(r0) -
+ sphere.GetRadius() * sphere.GetRadius();
+
+ auto const discriminant = alpha * alpha - 4 * beta * v0.squaredNorm();
+
+ //~ std::cout << "discriminant: " << discriminant << std::endl;
+ //~ std::cout << "alpha: " << alpha << std::endl;
+ //~ std::cout << "beta: " << beta << std::endl;
+
+ if (discriminant.magnitude() > 0) {
+ (-alpha - sqrt(discriminant)) / (2 * v0.squaredNorm());
+ return std::make_pair((-alpha - sqrt(discriminant)) / (2 * v0.squaredNorm()),
+ (-alpha + sqrt(discriminant)) / (2 * v0.squaredNorm()));
+ } else {
+ return {};
+ }
+ }
+
+ template <class Stack, class Trajectory>
+ TrackingLine<Stack, Trajectory>::TrackingLine(
+ corsika::environment::Environment const& pEnv)
+ : fEnvironment(pEnv) {}
+
+ template <class Stack, class Trajectory>
+ void TrackingLine<Stack, Trajectory>::Init() {}
+
+ template <class Stack, class Trajectory>
+ Trajectory TrackingLine<Stack, Trajectory>::GetTrack(Particle const& p) {
+ using std::cout;
+ using std::endl;
+ 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();
+ std::cout << "TrackingLine pid: " << p.GetPID() << " , E = " << p.GetEnergy() / 1_GeV
+ << " GeV" << std::endl;
+ std::cout << "TrackingLine pos: " << currentPosition.GetCoordinates() << std::endl;
+ std::cout << "TrackingLine E: " << p.GetEnergy() / 1_GeV << " GeV" << std::endl;
+ std::cout << "TrackingLine p: " << p.GetMomentum().GetComponents() / 1_GeV
+ << " GeV " << std::endl;
+ std::cout << "TrackingLine v: " << velocity.GetComponents() << std::endl;
+
+ // to do: include effect of magnetic field
+ geometry::Line line(currentPosition, velocity);
+
+ auto const* currentVolumeNode =
+ fEnvironment.GetUniverse()->GetContainingNode(currentPosition);
+ auto const& children = currentVolumeNode->GetChildNodes();
+ auto const& excluded = currentVolumeNode->GetExcludedNodes();
+
+ std::vector<TimeType> intersectionTimes;
+
+ auto addIfIntersects = [&](auto& 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;
+ if (t1.magnitude() >= 0)
+ intersectionTimes.push_back(t1);
+ else if (t2.magnitude() >= 0)
+ intersectionTimes.push_back(t2);
+ }
+ };
+
+ for (auto const& child : children) { addIfIntersects(*child); }
+
+ for (auto const* child : excluded) { addIfIntersects(*child); }
+
+ addIfIntersects(*currentVolumeNode);
+
+ auto const minIter =
+ std::min_element(intersectionTimes.cbegin(), intersectionTimes.cend());
+
+ TimeType min;
+
+ if (minIter == intersectionTimes.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;
+ }
+
+ std::cout << " t-intersect: " << min << std::endl;
+
+ return Trajectory(line, min);
+ }
+
+} // namespace corsika::process::tracking_line
+
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+using namespace corsika::setup;
+using Particle = Stack::ParticleType;
+// using Track = Trajectory;
+template class corsika::process::tracking_line::TrackingLine<setup::Stack,
+ setup::Trajectory>;
+
+#include "testTrackingLineStack.h"
+template class corsika::process::tracking_line::TrackingLine<DummyStack,
+ setup::Trajectory>;
diff --git a/Processes/TrackingLine/TrackingLine.h b/Processes/TrackingLine/TrackingLine.h
index 02126ee1..b5e715cd 100644
--- a/Processes/TrackingLine/TrackingLine.h
+++ b/Processes/TrackingLine/TrackingLine.h
@@ -12,29 +12,25 @@
#ifndef _include_corsika_processes_TrackingLine_h_
#define _include_corsika_processes_TrackingLine_h_
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/QuantityVector.h>
-#include <corsika/geometry/Sphere.h>
-#include <corsika/geometry/Vector.h>
-
#include <corsika/units/PhysicalUnits.h>
-
-#include <corsika/environment/Environment.h>
-#include <corsika/setup/SetupStack.h>
-#include <corsika/setup/SetupTrajectory.h>
-
-#include <algorithm>
-#include <iostream>
+#include <map> // for std::pair
#include <optional>
-#include <stdexcept>
-#include <utility>
+
+namespace corsika::environment {
+ class Environment;
+}
+namespace corsika::geometry {
+ class Line;
+ class Sphere;
+} // namespace corsika::geometry
namespace corsika::process {
namespace tracking_line {
- template <typename Stack>
+ template <typename Stack, typename Trajectory>
class TrackingLine { //
+
using Particle = typename Stack::ParticleType;
corsika::environment::Environment const& fEnvironment;
@@ -42,104 +38,13 @@ namespace corsika::process {
public:
std::optional<std::pair<corsika::units::si::TimeType, corsika::units::si::TimeType>>
TimeOfIntersection(corsika::geometry::Line const& line,
- geometry::Sphere const& sphere) {
- using namespace corsika::units::si;
- auto const& cs = fEnvironment.GetCoordinateSystem();
- geometry::Point const origin(cs, 0_m, 0_m, 0_m);
-
- auto const r0 = (line.GetR0() - origin);
- auto const v0 = line.GetV0();
- auto const c0 = (sphere.GetCenter() - origin);
-
- auto const alpha = r0.dot(v0) - 2 * v0.dot(c0);
- auto const beta = c0.squaredNorm() + r0.squaredNorm() + 2 * c0.dot(r0) -
- sphere.GetRadius() * sphere.GetRadius();
-
- auto const discriminant = alpha * alpha - 4 * beta * v0.squaredNorm();
-
- //~ std::cout << "discriminant: " << discriminant << std::endl;
- //~ std::cout << "alpha: " << alpha << std::endl;
- //~ std::cout << "beta: " << beta << std::endl;
-
- if (discriminant.magnitude() > 0) {
- (-alpha - sqrt(discriminant)) / (2 * v0.squaredNorm());
- return std::make_pair((-alpha - sqrt(discriminant)) / (2 * v0.squaredNorm()),
- (-alpha + sqrt(discriminant)) / (2 * v0.squaredNorm()));
- } else {
- return {};
- }
- }
-
- TrackingLine(corsika::environment::Environment const& pEnv)
- : fEnvironment(pEnv) {}
- void Init() {}
-
- auto GetTrack(Particle const& p) {
- using std::cout;
- using std::endl;
- 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();
- std::cout << "TrackingLine pid: " << p.GetPID()
- << " , E = " << p.GetEnergy() / 1_GeV << " GeV" << std::endl;
- std::cout << "TrackingLine pos: " << currentPosition.GetCoordinates()
- << std::endl;
- std::cout << "TrackingLine E: " << p.GetEnergy() / 1_GeV << " GeV" << std::endl;
- std::cout << "TrackingLine p: " << p.GetMomentum().GetComponents() / 1_GeV
- << " GeV " << std::endl;
- std::cout << "TrackingLine v: " << velocity.GetComponents() << std::endl;
-
- // to do: include effect of magnetic field
- geometry::Line line(currentPosition, velocity);
-
- auto const* currentVolumeNode =
- fEnvironment.GetUniverse()->GetContainingNode(currentPosition);
- auto const& children = currentVolumeNode->GetChildNodes();
- auto const& excluded = currentVolumeNode->GetExcludedNodes();
-
- std::vector<TimeType> intersectionTimes;
-
- auto addIfIntersects = [&](auto& 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;
- if (t1.magnitude() >= 0)
- intersectionTimes.push_back(t1);
- else if (t2.magnitude() >= 0)
- intersectionTimes.push_back(t2);
- }
- };
-
- for (auto const& child : children) { addIfIntersects(*child); }
-
- for (auto const* child : excluded) { addIfIntersects(*child); }
-
- addIfIntersects(*currentVolumeNode);
-
- auto const minIter =
- std::min_element(intersectionTimes.cbegin(), intersectionTimes.cend());
-
- TimeType min;
-
- if (minIter == intersectionTimes.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;
- }
-
- std::cout << " t-intersect: " << min << std::endl;
-
- return geometry::Trajectory<corsika::geometry::Line>(line, min);
- }
+ geometry::Sphere const& sphere);
+
+ TrackingLine(corsika::environment::Environment const& pEnv);
+
+ void Init();
+
+ Trajectory GetTrack(Particle const& p);
};
} // namespace tracking_line
diff --git a/Processes/TrackingLine/testTrackingLine.cc b/Processes/TrackingLine/testTrackingLine.cc
index 7dde7aec..00ad3e8e 100644
--- a/Processes/TrackingLine/testTrackingLine.cc
+++ b/Processes/TrackingLine/testTrackingLine.cc
@@ -9,11 +9,12 @@
* 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/process/tracking_line/TrackingLine.h>
-
#include <corsika/geometry/Point.h>
#include <corsika/geometry/Sphere.h>
#include <corsika/geometry/Vector.h>
@@ -34,33 +35,11 @@ using namespace corsika::geometry;
using namespace std;
using namespace corsika::units::si;
-typedef corsika::units::si::hepmomentum_d MOMENTUM;
-
-struct DummyParticle {
- HEPEnergyType fEnergy;
- Vector<MOMENTUM> fMomentum;
- Point fPosition;
-
- DummyParticle(HEPEnergyType pEnergy, Vector<MOMENTUM> pMomentum, Point pPosition)
- : fEnergy(pEnergy)
- , fMomentum(pMomentum)
- , fPosition(pPosition) {}
-
- auto GetEnergy() const { return fEnergy; }
- auto GetMomentum() const { return fMomentum; }
- auto GetPosition() const { return fPosition; }
- auto GetPID() const { return corsika::particles::Code::Unknown; }
-};
-
-struct DummyStack {
- using ParticleType = DummyParticle;
-};
-
TEST_CASE("TrackingLine") {
corsika::environment::Environment env; // dummy environment
auto const& cs = env.GetCoordinateSystem();
- tracking_line::TrackingLine<DummyStack> tracking(env);
+ tracking_line::TrackingLine<DummyStack, setup::Trajectory> tracking(env);
SECTION("intersection with sphere") {
Point const origin(cs, {0_m, 0_m, 0_m});
@@ -70,7 +49,7 @@ TEST_CASE("TrackingLine") {
0_m / second, 1_m / second);
Line line(origin, v);
- geometry::Trajectory<Line> traj(line, 12345_s);
+ setup::Trajectory traj(line, 12345_s);
auto const opt =
tracking.TimeOfIntersection(traj, Sphere(Point(cs, {0_m, 0_m, 10_m}), 1_m));
diff --git a/Processes/TrackingLine/testTrackingLineStack.cc b/Processes/TrackingLine/testTrackingLineStack.cc
new file mode 100644
index 00000000..0a718634
--- /dev/null
+++ b/Processes/TrackingLine/testTrackingLineStack.cc
@@ -0,0 +1,26 @@
+#ifndef _include_process_trackinling_teststack_h_
+#define _include_process_trackinling_teststack_h_
+
+typedef corsika::units::si::hepmomentum_d MOMENTUM;
+
+struct DummyParticle {
+ HEPEnergyType fEnergy;
+ Vector<MOMENTUM> fMomentum;
+ Point fPosition;
+
+ DummyParticle(HEPEnergyType pEnergy, Vector<MOMENTUM> pMomentum, Point pPosition)
+ : fEnergy(pEnergy)
+ , fMomentum(pMomentum)
+ , fPosition(pPosition) {}
+
+ auto GetEnergy() const { return fEnergy; }
+ auto GetMomentum() const { return fMomentum; }
+ auto GetPosition() const { return fPosition; }
+ auto GetPID() const { return corsika::particles::Code::Unknown; }
+};
+
+struct DummyStack {
+ using ParticleType = DummyParticle;
+};
+
+#endif
diff --git a/Processes/TrackingLine/testTrackingLineStack.h b/Processes/TrackingLine/testTrackingLineStack.h
new file mode 100644
index 00000000..6d8a46c1
--- /dev/null
+++ b/Processes/TrackingLine/testTrackingLineStack.h
@@ -0,0 +1,33 @@
+#ifndef _include_process_trackinling_teststack_h_
+#define _include_process_trackinling_teststack_h_
+
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/Vector.h>
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/units/PhysicalUnits.h>
+
+typedef corsika::units::si::hepmomentum_d MOMENTUM;
+
+struct DummyParticle {
+ corsika::units::si::HEPEnergyType fEnergy;
+ corsika::geometry::Vector<MOMENTUM> fMomentum;
+ corsika::geometry::Point fPosition;
+
+ DummyParticle(corsika::units::si::HEPEnergyType pEnergy,
+ corsika::geometry::Vector<MOMENTUM> pMomentum,
+ corsika::geometry::Point pPosition)
+ : fEnergy(pEnergy)
+ , fMomentum(pMomentum)
+ , fPosition(pPosition) {}
+
+ auto GetEnergy() const { return fEnergy; }
+ auto GetMomentum() const { return fMomentum; }
+ auto GetPosition() const { return fPosition; }
+ auto GetPID() const { return corsika::particles::Code::Unknown; }
+};
+
+struct DummyStack {
+ using ParticleType = DummyParticle;
+};
+
+#endif
diff --git a/Setup/SetupStack.h b/Setup/SetupStack.h
index 397f7a35..9f42d76c 100644
--- a/Setup/SetupStack.h
+++ b/Setup/SetupStack.h
@@ -68,13 +68,13 @@ public:
void SetParticleData(const std::tuple<const corsika::environment::BaseNodeType*> v) {
SetNode(std::get<0>(v));
}
- void SetParticleData(
- GeometryDataInterface<T>& parent,
- const std::tuple<const corsika::environment::BaseNodeType*>) {
+ void SetParticleData(GeometryDataInterface<T>& parent,
+ const std::tuple<const corsika::environment::BaseNodeType*>) {
SetNode(parent.GetNode()); // copy Node from parent particle!
}
void SetParticleData() { SetNode(nullptr); }
- void SetParticleData(GeometryDataInterface<T>& parent) { SetNode(parent.GetNode()); // copy Node from parent particle!
+ void SetParticleData(GeometryDataInterface<T>& parent) {
+ SetNode(parent.GetNode()); // copy Node from parent particle!
}
void SetNode(const corsika::environment::BaseNodeType* v) {
GetStackData().SetNode(GetIndex(), v);
--
GitLab