IAP GITLAB

Commit dbf6e995 authored by Ralf Ulrich's avatar Ralf Ulrich

made leap-frog more stabel and robust, also in extrem magnetic field cases

parent eac7f6c0
Pipeline #2890 failed with stages
in 0 seconds
......@@ -8,7 +8,6 @@
#include <corsika/cascade/Cascade.h>
#include <corsika/process/ProcessSequence.h>
#include <corsika/process/tracking_line/Tracking.h>
#include <corsika/setup/SetupEnvironment.h>
#include <corsika/setup/SetupStack.h>
......@@ -123,7 +122,7 @@ int main() {
universe.AddChild(std::move(world));
// setup processes, decays and interactions
tracking_line::Tracking tracking;
setup::Tracking tracking;
random::RNGManager::GetInstance().RegisterRandomStream("sibyll");
process::sibyll::Interaction sibyll;
......
......@@ -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/Tracking.h>
#include <corsika/setup/SetupEnvironment.h>
#include <corsika/setup/SetupStack.h>
......@@ -135,7 +134,7 @@ int main() {
}
// setup processes, decays and interactions
tracking_line::Tracking tracking;
setup::Tracking tracking;
stack_inspector::StackInspector<setup::Stack> stackInspect(1, true, E0);
random::RNGManager::GetInstance().RegisterRandomStream("sibyll");
......
......@@ -10,7 +10,6 @@
#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/Tracking.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
......@@ -121,7 +120,7 @@ int main() {
}
// setup processes, decays and interactions
tracking_line::Tracking tracking;
setup::Tracking tracking;
stack_inspector::StackInspector<setup::Stack> stackInspect(1, true, E0);
random::RNGManager::GetInstance().RegisterRandomStream("sibyll");
......
......@@ -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/Tracking.h>
#include <corsika/random/RNGManager.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
......@@ -157,7 +156,7 @@ int main(int argc, char** argv) {
auto sequence = process::sequence(proposalCounted, em_continuous, longprof, cut,
observationLevel, trackWriter);
// define air shower object, run simulation
tracking_line::Tracking tracking;
setup::Tracking tracking;
cascade::Cascade EAS(env, tracking, sequence, stack);
// to fix the point of first interaction, uncomment the following two lines:
......
......@@ -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/Tracking.h>
#include <corsika/process/urqmd/UrQMD.h>
#include <corsika/random/RNGManager.h>
#include <corsika/setup/SetupStack.h>
......@@ -242,7 +241,7 @@ int main(int argc, char** argv) {
eLoss, cut, conex, longprof, observationLevel);
// define air shower object, run simulation
tracking_line::Tracking tracking;
setup::Tracking tracking;
cascade::Cascade EAS(env, tracking, sequence, stack);
// to fix the point of first interaction, uncomment the following two lines:
......
......@@ -38,9 +38,9 @@
#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/Tracking.h>
#include <corsika/process/urqmd/UrQMD.h>
#include <corsika/random/RNGManager.h>
#include <corsika/setup/SetupStack.h>
......@@ -85,7 +85,7 @@ using MyExtraEnv = environment::MediumPropertyModel<environment::UniformMagnetic
int main(int argc, char** argv) {
logging::SetLevel(logging::level::info);
logging::SetLevel(logging::level::trace);
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, 50_uT, 0_T});
geometry::Vector{rootCS, 0_T, 5000_mT, 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
......@@ -255,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);
stack_inspector::StackInspector<setup::Stack> stackInspect(1000, false, E0);
auto sequence =
process::sequence(hadronSequence, reset_particle_mass, decaySequence,
proposalCounted, em_continuous, cut, observationLevel, longprof);
process::sequence(stackInspect, hadronSequence, reset_particle_mass, decaySequence,
proposalCounted, em_continuous, cut, trackWriter, observationLevel, longprof);
// define air shower object, run simulation
tracking_line::Tracking tracking;
setup::Tracking tracking;
cascade::Cascade EAS(env, tracking, sequence, stack);
// to fix the point of first interaction, uncomment the following two lines:
......
......@@ -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,7 +47,7 @@ namespace corsika::environment {
*/
// clang-format on
units::si::GrammageType IntegratedGrammage(
geometry::LineTrajectory 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(); }
......@@ -80,7 +80,7 @@ namespace corsika::environment {
*/
// clang-format on
units::si::LengthType ArclengthFromGrammage(
geometry::LineTrajectory 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.GetDirection(0).dot(vAxis).magnitude();
auto const rhoStart = GetImplementation().GetMassDensity(vLine.GetLine().GetR0());
......
......@@ -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 {
......
......@@ -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,13 +52,13 @@ namespace corsika::environment {
NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
units::si::GrammageType IntegratedGrammage(geometry::LineTrajectory const& vLine,
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::LineTrajectory const& vLine,
setup::Trajectory const& vLine,
units::si::GrammageType vGrammage) const override {
return Base::ArclengthFromGrammage(vLine, vGrammage, fAxis);
}
......
......@@ -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::LineTrajectory 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::LineTrajectory const&,
corsika::setup::Trajectory const&,
corsika::units::si::GrammageType pGrammage) const override {
return pGrammage / fDensity;
}
......
......@@ -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::LineTrajectory const&,
corsika::setup::Trajectory const&,
corsika::units::si::LengthType) const = 0;
virtual corsika::units::si::LengthType ArclengthFromGrammage(
corsika::geometry::LineTrajectory const&,
corsika::setup::Trajectory const&,
corsika::units::si::GrammageType) const = 0;
virtual NuclearComposition const& GetNuclearComposition() const = 0;
......
......@@ -41,13 +41,13 @@ namespace corsika::environment {
NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
corsika::units::si::GrammageType IntegratedGrammage(
corsika::geometry::LineTrajectory 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::LineTrajectory const& pLine,
corsika::setup::Trajectory const& pLine,
corsika::units::si::GrammageType pGrammage) const override {
return fDensityFunction.ArclengthFromGrammage(pLine, pGrammage);
}
......
......@@ -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,7 +20,7 @@ namespace corsika::environment {
auto const& GetImplementation() const { return *static_cast<TDerived const*>(this); }
public:
auto IntegrateGrammage(corsika::geometry::LineTrajectory const& line,
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),
......@@ -28,7 +28,7 @@ namespace corsika::environment {
return (c0 + 0.5 * c1 * length) * length;
}
auto ArclengthFromGrammage(corsika::geometry::LineTrajectory const& line,
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),
......@@ -37,7 +37,7 @@ namespace corsika::environment {
return (1 - 0.5 * grammage * c1 / (c0 * c0)) * grammage / c0;
}
auto MaximumLength(corsika::geometry::LineTrajectory 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(
......
......@@ -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_; }
......
......@@ -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_;
......
......@@ -55,14 +55,14 @@ namespace corsika::environment {
NuclearComposition const& GetNuclearComposition() const override { return nuclComp_; }
units::si::GrammageType IntegratedGrammage(
geometry::LineTrajectory const& line,
setup::Trajectory const& line,
units::si::LengthType l) const override {
auto const axis = (line.GetLine().GetR0() - Base::fP0).normalized();
return Base::IntegratedGrammage(line, l, axis);
}
units::si::LengthType ArclengthFromGrammage(
geometry::LineTrajectory const& line,
setup::Trajectory const& line,
units::si::GrammageType grammage) const override {
auto const axis = (line.GetLine().GetR0() - Base::fP0).normalized();
return Base::ArclengthFromGrammage(line, grammage, axis);
......
......@@ -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}));
LineTrajectory 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}));
LineTrajectory 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,7 +82,7 @@ TEST_CASE("FlatExponential") {
SECTION("escape grammage") {
Line const line(gOrigin, Vector<SpeedType::dimension_type>(
gCS, {0_m / second, 0_m / second, -5_m / second}));
LineTrajectory const trajectory(line, tEnd);
setup::Trajectory const trajectory = setup::testing::make_track<setup::Trajectory>(line, tEnd);
GrammageType const escapeGrammage = rho0 * lambda;
......@@ -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}));
LineTrajectory 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}));
LineTrajectory 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;
LineTrajectory 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
LineTrajectory 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
LineTrajectory 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
LineTrajectory 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));
......
......@@ -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);
......
......@@ -111,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());
......@@ -202,21 +202,20 @@ namespace corsika::cascade {
next_interact);
// determine the maximum geometric step length
LengthType const distance_max = process_sequence_.MaxStepLength(vParticle, step);
C8LOG_DEBUG("distance_max={} m", distance_max / 1_m);
LengthType const continuous_max_dist = process_sequence_.MaxStepLength(vParticle, step);
// take minimum of geometry, interaction, decay for next step
auto min_distance =
std::min({distance_interact, distance_decay, distance_max, geomMaxLength});
std::min({distance_interact, distance_decay, continuous_max_dist, geomMaxLength});
C8LOG_DEBUG(
"transport particle by : {} m "
"Medium transition after: {} m "
"Decay after: {} m "
"Interaction after: {} m"
"Continuous limit: {} m",
"Interaction after: {} m "
"Continuous limit: {} m ",
min_distance / 1_m, geomMaxLength / 1_m, distance_decay / 1_m,
distance_interact / 1_m, distance_max / 1_m);
distance_interact / 1_m, continuous_max_dist / 1_m);
// here the particle is actually moved along the trajectory to new position:
step.SetLength(min_distance);
......@@ -248,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
......@@ -261,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 &&
......@@ -289,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);
}
}
}
......
......@@ -47,7 +47,9 @@ namespace corsika::geometry {
/**
* \param theLine The geometric \sa Line object that represents a straight-line
* connection \param timeLength The time duration to traverse the straight trajectory
* 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)
......@@ -59,9 +61,15 @@ namespace corsika::geometry {
/**
* \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
* 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
*/