From fd16e22ae23c1122cdf0d5cb92bb0251afb5647d Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Tue, 24 Nov 2020 06:18:30 +0100
Subject: [PATCH] made many small details more robust, made a lot of testing
---
Documentation/Examples/boundary_example.cc | 2 +-
Documentation/Examples/cascade_example.cc | 6 +-
.../Examples/cascade_proton_example.cc | 25 ++--
Documentation/Examples/em_shower.cc | 4 +-
Documentation/Examples/hybrid_MC.cc | 2 +-
Documentation/Examples/vertical_EAS.cc | 10 +-
Processes/TrackingLeapFrogCurved/Tracking.h | 109 +++++++++++++-----
Processes/TrackingLeapFrogStraight/Tracking.h | 89 +++++++++-----
Processes/TrackingLine/Tracking.h | 2 +
Setup/SetupTrajectory.h | 31 ++++-
10 files changed, 196 insertions(+), 84 deletions(-)
diff --git a/Documentation/Examples/boundary_example.cc b/Documentation/Examples/boundary_example.cc
index c40823c14..ed9d7bd43 100644
--- a/Documentation/Examples/boundary_example.cc
+++ b/Documentation/Examples/boundary_example.cc
@@ -101,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<
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 755c7c7c1..8b9b5b8ad 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -57,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;
@@ -72,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<
@@ -110,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) {
diff --git a/Documentation/Examples/cascade_proton_example.cc b/Documentation/Examples/cascade_proton_example.cc
index 909cad004..b36239219 100644
--- a/Documentation/Examples/cascade_proton_example.cc
+++ b/Documentation/Examples/cascade_proton_example.cc
@@ -10,6 +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/energy_loss/EnergyLoss.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
@@ -17,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>
@@ -59,7 +61,7 @@ using namespace corsika::units::si;
//
int main() {
- logging::SetLevel(logging::level::trace);
+ logging::SetLevel(logging::level::info);
std::cout << "cascade_proton_example" << std::endl;
@@ -73,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;
@@ -97,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);
@@ -112,11 +115,10 @@ 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
@@ -130,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 31298aa28..07c4cfb7b 100644
--- a/Documentation/Examples/em_shower.cc
+++ b/Documentation/Examples/em_shower.cc
@@ -71,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
diff --git a/Documentation/Examples/hybrid_MC.cc b/Documentation/Examples/hybrid_MC.cc
index 75ed47a19..2216e4b33 100644
--- a/Documentation/Examples/hybrid_MC.cc
+++ b/Documentation/Examples/hybrid_MC.cc
@@ -101,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
diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc
index 834989466..6fbed4a73 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -85,7 +85,7 @@ using MyExtraEnv = environment::MediumPropertyModel<environment::UniformMagnetic
int main(int argc, char** argv) {
- logging::SetLevel(logging::level::trace);
+ logging::SetLevel(logging::level::info);
C8LOG_INFO("vertical_EAS");
@@ -110,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, 5000_mT, 0_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
@@ -256,9 +256,9 @@ int main(int argc, char** argv) {
EnergySwitch(55_GeV));
auto decaySequence = process::sequence(decayPythia, decaySibyll);
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);
+ auto sequence = process::sequence(stackInspect, hadronSequence, reset_particle_mass,
+ decaySequence, proposalCounted, em_continuous, cut,
+ trackWriter, observationLevel, longprof);
// define air shower object, run simulation
setup::Tracking tracking;
diff --git a/Processes/TrackingLeapFrogCurved/Tracking.h b/Processes/TrackingLeapFrogCurved/Tracking.h
index 9993f7b70..da636d656 100644
--- a/Processes/TrackingLeapFrogCurved/Tracking.h
+++ b/Processes/TrackingLeapFrogCurved/Tracking.h
@@ -39,6 +39,8 @@ namespace corsika::process {
* \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,
@@ -50,7 +52,7 @@ namespace corsika::process {
} // charge of the particle
const int chargeNumber = particle.GetChargeNumber();
auto const* currentLogicalVolumeNode = particle.GetNode();
- auto magneticfield =
+ MagneticFieldVector const& magneticfield =
currentLogicalVolumeNode->GetModelProperties().GetMagneticField(
particle.GetPosition());
geometry::Vector<SpeedType::dimension_type> velocity =
@@ -85,6 +87,9 @@ namespace corsika::process {
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;
@@ -111,23 +116,34 @@ namespace corsika::process {
// maximum step-length since we need to follow curved
// trajectories segment-wise -- at least if we don't employ concepts as "Helix
// Trajectories" or similar
- const auto& magneticfield =
+ MagneticFieldVector const& magneticfield =
volumeNode.GetModelProperties().GetMagneticField(position);
- const auto magnitudeB = magneticfield.norm();
- const int chargeNumber = particle.GetChargeNumber();
- auto const momentumVerticalMag =
- particle.GetMomentum() -
- particle.GetMomentum().parallelProjectionOnto(magneticfield);
+ 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 =
- (chargeNumber == 0 || magnitudeB == 0_T
- ? std::numeric_limits<TimeType::value_type>::infinity() * 1_m
- : momentumVerticalMag.norm() * 1_V /
- (corsika::units::constants::c * abs(chargeNumber) * magnitudeB *
- 1_eV));
+ (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: {} m = {} s", gyroradius, steplimit,
+ C8LOG_DEBUG("gyroradius {}, steplimit: {} = {}", gyroradius, steplimit,
steplimit_time);
// traverse the environment volume tree and find next
@@ -148,39 +164,48 @@ namespace corsika::process {
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();
- const auto& magneticfield = medium.GetMagneticField(position);
+ 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
- if (chargeNumber == 0 || magneticfield.norm() == 0_T) {
+ 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 geometry::Vector<SpeedType::dimension_type> velocity =
- particle.GetMomentum() / particle.GetEnergy() * corsika::units::constants::c;
const auto absVelocity = velocity.norm();
auto energy = particle.GetEnergy();
auto k = chargeNumber * corsika::units::constants::cSquared * 1_eV /
(absVelocity * energy * 1_V);
- const geometry::Vector<dimensionless_d> directionBefore =
- velocity.normalized(); // determine steplength to next volume
+ 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 * (directionBefore.cross(magneticfield)).GetSquaredNorm() * k * k);
+ 4 / (1_m * 1_m * denom);
const double b = directionBefore.dot(position - sphere.GetCenter()) * 8 /
- ((directionBefore.cross(magneticfield)).GetSquaredNorm() * k *
- k * 1_m * 1_m * 1_m);
+ (denom * 1_m * 1_m * 1_m);
const double c = ((position - sphere.GetCenter()).GetSquaredNorm() -
(sphere.GetRadius() * sphere.GetRadius())) *
- 4 /
- ((directionBefore.cross(magneticfield)).GetSquaredNorm() * k *
- k * 1_m * 1_m * 1_m * 1_m);
+ 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;
@@ -252,6 +277,38 @@ namespace corsika::process {
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
diff --git a/Processes/TrackingLeapFrogStraight/Tracking.h b/Processes/TrackingLeapFrogStraight/Tracking.h
index b88a6fc46..efa46981f 100644
--- a/Processes/TrackingLeapFrogStraight/Tracking.h
+++ b/Processes/TrackingLeapFrogStraight/Tracking.h
@@ -52,6 +52,19 @@ namespace corsika::process {
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;
@@ -72,39 +85,9 @@ namespace corsika::process {
typedef decltype(particle.GetNode()) node_type;
const node_type volumeNode = particle.GetNode();
- auto magneticfield =
- volumeNode->GetModelProperties().GetMagneticField(initialPosition);
-
- // charge of the particle
- const int chargeNumber = particle.GetChargeNumber();
- const auto magnitudeB = magneticfield.GetNorm();
- C8LOG_DEBUG("field={} uT, chargeNumber={}, magnitudeB={} uT",
- magneticfield.GetComponents() / 1_uT, chargeNumber, magnitudeB / 1_T);
-
- // we need to limit maximum step-length since we absolutely
- // need to follow strongly curved trajectories segment-wise,
- // at least if we don't employ concepts as "Helix
- // Trajectories" or similar
- auto const momentumVerticalMag =
- particle.GetMomentum() -
- particle.GetMomentum().parallelProjectionOnto(magneticfield);
- bool const no_deflection = chargeNumber == 0 || magnitudeB == 0_T;
- LengthType const gyroradius =
- (no_deflection ? std::numeric_limits<TimeType::value_type>::infinity() * 1_m
- : momentumVerticalMag.norm() * 1_V /
- (corsika::units::constants::c * abs(chargeNumber) *
- magnitudeB * 1_eV));
- const double maxRadians = 0.01;
- const LengthType steplimit = 2 * cos(maxRadians) * sin(maxRadians) * gyroradius;
- C8LOG_DEBUG("gyroradius {}, Steplimit: {}", gyroradius, steplimit);
-
- // calculate first halve step for "steplimit"
- const auto initialMomentum = particle.GetMomentum();
- const auto absMomentum = initialMomentum.norm();
- const auto absVelocity = initialVelocity.norm();
- const geometry::Vector<dimensionless_d> direction = initialVelocity.normalized();
// check if particle is moving at all
+ const auto absVelocity = initialVelocity.norm();
if (absVelocity * 1_s == 0_m) {
return std::make_tuple(
geometry::LineTrajectory(geometry::Line(initialPosition, initialVelocity),
@@ -112,6 +95,15 @@ namespace corsika::process {
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);
@@ -128,9 +120,39 @@ namespace corsika::process {
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) / 2;
+ std::min(steplimit, initialTrackLength * firstFraction_);
C8LOG_DEBUG("first halve step length {}, steplimit={}, initialTrackLength={}",
firstHalveSteplength, steplimit, initialTrackLength);
@@ -203,6 +225,9 @@ namespace corsika::process {
new_direction_normalized * absVelocity), // trajectory
(switch_volume ? finalTrackNextVolume : volumeNode));
}
+
+ protected:
+ double firstFraction_;
};
} // namespace tracking_leapfrog_straight
diff --git a/Processes/TrackingLine/Tracking.h b/Processes/TrackingLine/Tracking.h
index d9690ec0d..0dfac370a 100644
--- a/Processes/TrackingLine/Tracking.h
+++ b/Processes/TrackingLine/Tracking.h
@@ -35,6 +35,7 @@ namespace corsika::process {
class Tracking : public corsika::process::tracking::Intersect<Tracking> {
public:
+
template <typename TParticle>
auto GetTrack(TParticle const& particle) {
using namespace corsika::units::si;
@@ -115,6 +116,7 @@ namespace corsika::process {
1_s
: n.dot(delta) / c);
}
+
};
} // namespace tracking_line
diff --git a/Setup/SetupTrajectory.h b/Setup/SetupTrajectory.h
index a8b86e030..2a42db7a6 100644
--- a/Setup/SetupTrajectory.h
+++ b/Setup/SetupTrajectory.h
@@ -20,14 +20,39 @@
namespace corsika::setup {
- // Note: Tracking and Trajectory must fit together !
+ /**
+ Note/Warning: Tracking and Trajectory must fit together !
+
+ 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.
+
+ 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.
+
+ tracking_line::Tracking is a pure straight tracker. It is based on
+ LineTrajectory.
+ */
typedef corsika::process::tracking_leapfrog_curved::Tracking Tracking;
- // tracking_leapfrog_straight::Tracking tracking;
+ //typedef corsika::process::tracking_leapfrog_straight::Tracking Tracking;
+ //typedef corsika::process::tracking_line::Tracking Tracking;
/// definition of Trajectory base class, to be used in tracking and cascades
- // typedef corsika::geometry::LineTrajectory Trajectory;
+ //typedef corsika::geometry::LineTrajectory Trajectory;
typedef corsika::geometry::LeapFrogTrajectory Trajectory;
+ /**
+
+ The following section is for unit testing only. Eventually it should
+ be moved to "tests".
+
+
+ */
+
namespace testing {
template <typename TTrack>
--
GitLab