diff --git a/corsika/detail/framework/geometry/Cubic.inl b/corsika/detail/framework/geometry/Cubic.inl
index be944049e043aee323e00a2e604e0aa16114081f..5e64f3abaf8fdbf8cc88ad5b9bfc75529312b019 100644
--- a/corsika/detail/framework/geometry/Cubic.inl
+++ b/corsika/detail/framework/geometry/Cubic.inl
@@ -1,20 +1,19 @@
namespace corsika {
-inline bool Cubic::contains(Point const &p) const {
- if ((abs(p.getX(cs_)) < x_) && (abs(p.getY(cs_)) < y_) &&
- (abs(p.getZ(cs_)) < z_))
- return true;
- else
- return false;
-}
+ inline bool Cubic::contains(Point const& p) const {
+ if ((abs(p.getX(cs_)) < x_) && (abs(p.getY(cs_)) < y_) && (abs(p.getZ(cs_)) < z_))
+ return true;
+ else
+ return false;
+ }
-inline std::string Cubic::asString() const {
- std::ostringstream txt;
- txt << "center=" << center_ << ", x-axis=" << DirectionVector{cs_, {1, 0, 0}}
- << ", y-axis: " << DirectionVector{cs_, {0, 1, 0}}
- << ", z-axis: " << DirectionVector{cs_, {0, 0, 1}};
- return txt.str();
-}
+ inline std::string Cubic::asString() const {
+ std::ostringstream txt;
+ txt << "center=" << center_ << ", x-axis=" << DirectionVector{cs_, {1, 0, 0}}
+ << ", y-axis: " << DirectionVector{cs_, {0, 1, 0}}
+ << ", z-axis: " << DirectionVector{cs_, {0, 0, 1}};
+ return txt.str();
+ }
} // namespace corsika
\ No newline at end of file
diff --git a/corsika/detail/modules/ObservationCubic.inl b/corsika/detail/modules/ObservationCubic.inl
index 1967e3dd66375babfabb530d179a09f1c98cee64..1504c67bc454336bd8af9bb55d09ce4e784e0d16 100644
--- a/corsika/detail/modules/ObservationCubic.inl
+++ b/corsika/detail/modules/ObservationCubic.inl
@@ -8,142 +8,144 @@
namespace corsika {
-template <typename TTracking, typename TOutput>
-ObservationCubic<TTracking, TOutput>::ObservationCubic(
- Point const ¢er, CoordinateSystemPtr cs, LengthType const x,
- LengthType const y, LengthType const z, bool deleteOnHit)
- : Cubic(center, cs, x, y, z), deleteOnHit_(deleteOnHit), energy_(0_GeV),
- count_(0) {}
-
-template <typename TTracking, typename TOutput>
-template <typename TParticle, typename TTrajectory>
-inline ProcessReturn ObservationCubic<TTracking, TOutput>::doContinuous(
- TParticle &particle, TTrajectory & step, bool const stepLimit) {
- /*
- The current step did not yet reach the ObservationCubic, do nothing now and
- wait:
- */
- if (!stepLimit) {
- // @todo this is actually needed to fix small instabilities of the leap-frog
- // tracking: Note, this is NOT a general solution and should be clearly
- // revised with a more robust tracking. #ifdef DEBUG
- if (deleteOnHit_) {
- // since this is basically a bug, it cannot be tested LCOV_EXCL_START
- LengthType const check = 1_m; // TODO, do I need to check?
- if (check < 0_m) {
- CORSIKA_LOG_WARN("PARTICLE AVOIDED ObservationCubic {}", check);
- CORSIKA_LOG_WARN("Temporary fix: write and remove particle.");
+ template <typename TTracking, typename TOutput>
+ ObservationCubic<TTracking, TOutput>::ObservationCubic(
+ Point const& center, CoordinateSystemPtr cs, LengthType const x, LengthType const y,
+ LengthType const z, bool deleteOnHit)
+ : Cubic(center, cs, x, y, z)
+ , deleteOnHit_(deleteOnHit)
+ , energy_(0_GeV)
+ , count_(0) {}
+
+ template <typename TTracking, typename TOutput>
+ template <typename TParticle, typename TTrajectory>
+ inline ProcessReturn ObservationCubic<TTracking, TOutput>::doContinuous(
+ TParticle& particle, TTrajectory& step, bool const stepLimit) {
+ /*
+ The current step did not yet reach the ObservationCubic, do nothing now and
+ wait:
+ */
+ if (!stepLimit) {
+ // @todo this is actually needed to fix small instabilities of the leap-frog
+ // tracking: Note, this is NOT a general solution and should be clearly
+ // revised with a more robust tracking. #ifdef DEBUG
+ if (deleteOnHit_) {
+ // since this is basically a bug, it cannot be tested LCOV_EXCL_START
+ LengthType const check = 1_m; // TODO, do I need to check?
+ if (check < 0_m) {
+ CORSIKA_LOG_WARN("PARTICLE AVOIDED ObservationCubic {}", check);
+ CORSIKA_LOG_WARN("Temporary fix: write and remove particle.");
+ } else
+ return ProcessReturn::Ok;
+ // LCOV_EXCL_STOP
} else
+ // #endif
return ProcessReturn::Ok;
- // LCOV_EXCL_STOP
- } else
- // #endif
+ }
+
+ HEPEnergyType const energy = particle.getEnergy();
+ Point const pointOfIntersection = step.getPosition(1);
+ DirectionVector const dirction = particle.getDirection();
+
+ // add our particles to the output file stream
+ this->write(particle.getPID(), energy, pointOfIntersection.getX(cs_),
+ pointOfIntersection.getY(cs_), pointOfIntersection.getZ(cs_),
+ dirction.getX(cs_), dirction.getY(cs_), dirction.getZ(cs_),
+ particle.getTime());
+
+ CORSIKA_LOG_TRACE("Particle detected absorbed={}", deleteOnHit_);
+
+ if (deleteOnHit_) {
+ count_++;
+ energy_ += energy;
+ return ProcessReturn::ParticleAbsorbed;
+ } else {
return ProcessReturn::Ok;
+ }
}
- HEPEnergyType const energy = particle.getEnergy();
- Point const pointOfIntersection = step.getPosition(1);
- DirectionVector const dirction = particle.getDirection();
+ template <typename TTracking, typename TOutput>
+ template <typename TParticle, typename TTrajectory>
+ inline LengthType ObservationCubic<TTracking, TOutput>::getMaxStepLength(
+ TParticle const& particle, TTrajectory const& trajectory) {
+
+ CORSIKA_LOG_TRACE("getMaxStepLength, particle={}, pos={}, dir={}, cubic={}",
+ particle.asString(), particle.getPosition(),
+ particle.getDirection(), asString());
+
+ auto const intersection =
+ TTracking::intersect(particle, static_cast<Cubic const>(*this));
+
+ TimeType const timeOfIntersection = intersection.getEntry();
+ CORSIKA_LOG_TRACE("timeOfIntersection={}", timeOfIntersection);
+ if (timeOfIntersection < TimeType::zero()) {
+ return std::numeric_limits<double>::infinity() * 1_m;
+ }
+ if (timeOfIntersection > trajectory.getDuration()) {
+ return std::numeric_limits<double>::infinity() * 1_m;
+ }
+ double const fractionOfIntersection = timeOfIntersection / trajectory.getDuration();
+ CORSIKA_LOG_TRACE("ObservationCubic: getMaxStepLength dist={} m, pos={}",
+ trajectory.getLength(fractionOfIntersection) / 1_m,
+ trajectory.getPosition(fractionOfIntersection));
+ return trajectory.getLength(fractionOfIntersection);
+ }
- // add our particles to the output file stream
- this->write(particle.getPID(), energy, pointOfIntersection.getX(cs_),
- pointOfIntersection.getY(cs_), pointOfIntersection.getZ(cs_),
- dirction.getX(cs_), dirction.getY(cs_), dirction.getZ(cs_),
- particle.getTime());
+ template <typename TTracking, typename TOutput>
+ inline void ObservationCubic<TTracking, TOutput>::showResults() const {
+ CORSIKA_LOG_INFO(
+ " ******************************\n"
+ " ObservationCubic: \n"
+ " energy at cubic (GeV) : {}\n"
+ " no. of particles at cubic : {}\n"
+ " ******************************",
+ energy_ / 1_GeV, count_);
+ }
- CORSIKA_LOG_TRACE("Particle detected absorbed={}", deleteOnHit_);
+ template <typename TTracking, typename TOutput>
+ inline YAML::Node ObservationCubic<TTracking, TOutput>::getConfig() const {
+ using namespace units::si;
- if (deleteOnHit_) {
- count_++;
- energy_ += energy;
- return ProcessReturn::ParticleAbsorbed;
- } else {
- return ProcessReturn::Ok;
- }
-}
+ // construct the top-level node
+ YAML::Node node;
+
+ // basic info
+ node["type"] = "ObservationCubic";
+ node["units"] = "m"; // add default units for values
-template <typename TTracking, typename TOutput>
-template <typename TParticle, typename TTrajectory>
-inline LengthType ObservationCubic<TTracking, TOutput>::getMaxStepLength(
- TParticle const &particle, TTrajectory const &trajectory) {
+ // save each component in its native coordinate system
+ auto const root_cs = get_root_CoordinateSystem();
+ node["center"].push_back(center_.getX(root_cs) / 1_m);
+ node["center"].push_back(center_.getY(root_cs) / 1_m);
+ node["center"].push_back(center_.getZ(root_cs) / 1_m);
- CORSIKA_LOG_TRACE("getMaxStepLength, particle={}, pos={}, dir={}, cubic={}",
- particle.asString(), particle.getPosition(),
- particle.getDirection(), asString());
+ // the x-axis vector
+ DirectionVector const x_axis = DirectionVector{cs_, {1, 0, 0}};
+ node["x-axis"].push_back(x_axis.getX(root_cs).magnitude());
+ node["x-axis"].push_back(x_axis.getY(root_cs).magnitude());
+ node["x-axis"].push_back(x_axis.getZ(root_cs).magnitude());
- auto const intersection =
- TTracking::intersect(particle, static_cast<Cubic const>(*this));
+ // the y-axis vector
+ DirectionVector const y_axis = DirectionVector{cs_, {0, 1, 0}};
+ node["y-axis"].push_back(y_axis.getX(root_cs).magnitude());
+ node["y-axis"].push_back(y_axis.getY(root_cs).magnitude());
+ node["y-axis"].push_back(y_axis.getZ(root_cs).magnitude());
- TimeType const timeOfIntersection = intersection.getEntry();
- CORSIKA_LOG_TRACE("timeOfIntersection={}", timeOfIntersection);
- if (timeOfIntersection < TimeType::zero()) {
- return std::numeric_limits<double>::infinity() * 1_m;
+ // the x-axis vector
+ DirectionVector const z_axis = DirectionVector{cs_, {0, 0, 1}};
+ node["z-axis"].push_back(z_axis.getX(root_cs).magnitude());
+ node["z-axis"].push_back(z_axis.getY(root_cs).magnitude());
+ node["z-axis"].push_back(z_axis.getZ(root_cs).magnitude());
+
+ node["delete_on_hit"] = deleteOnHit_;
+
+ return node;
}
- if (timeOfIntersection > trajectory.getDuration()) {
- return std::numeric_limits<double>::infinity() * 1_m;
+
+ template <typename TTracking, typename TOutput>
+ inline void ObservationCubic<TTracking, TOutput>::reset() {
+ energy_ = 0_GeV;
+ count_ = 0;
}
- double const fractionOfIntersection =
- timeOfIntersection / trajectory.getDuration();
- CORSIKA_LOG_TRACE("ObservationCubic: getMaxStepLength dist={} m, pos={}",
- trajectory.getLength(fractionOfIntersection) / 1_m,
- trajectory.getPosition(fractionOfIntersection));
- return trajectory.getLength(fractionOfIntersection);
-}
-
-template <typename TTracking, typename TOutput>
-inline void ObservationCubic<TTracking, TOutput>::showResults() const {
- CORSIKA_LOG_INFO(" ******************************\n"
- " ObservationCubic: \n"
- " energy at cubic (GeV) : {}\n"
- " no. of particles at cubic : {}\n"
- " ******************************",
- energy_ / 1_GeV, count_);
-}
-
-template <typename TTracking, typename TOutput>
-inline YAML::Node ObservationCubic<TTracking, TOutput>::getConfig() const {
- using namespace units::si;
-
- // construct the top-level node
- YAML::Node node;
-
- // basic info
- node["type"] = "ObservationCubic";
- node["units"] = "m"; // add default units for values
-
- // save each component in its native coordinate system
- auto const root_cs = get_root_CoordinateSystem();
- node["center"].push_back(center_.getX(root_cs) / 1_m);
- node["center"].push_back(center_.getY(root_cs) / 1_m);
- node["center"].push_back(center_.getZ(root_cs) / 1_m);
-
- // the x-axis vector
- DirectionVector const x_axis = DirectionVector{cs_, {1, 0, 0}};
- node["x-axis"].push_back(x_axis.getX(root_cs).magnitude());
- node["x-axis"].push_back(x_axis.getY(root_cs).magnitude());
- node["x-axis"].push_back(x_axis.getZ(root_cs).magnitude());
-
- // the y-axis vector
- DirectionVector const y_axis = DirectionVector{cs_, {0, 1, 0}};
- node["y-axis"].push_back(y_axis.getX(root_cs).magnitude());
- node["y-axis"].push_back(y_axis.getY(root_cs).magnitude());
- node["y-axis"].push_back(y_axis.getZ(root_cs).magnitude());
-
- // the x-axis vector
- DirectionVector const z_axis = DirectionVector{cs_, {0, 0, 1}};
- node["z-axis"].push_back(z_axis.getX(root_cs).magnitude());
- node["z-axis"].push_back(z_axis.getY(root_cs).magnitude());
- node["z-axis"].push_back(z_axis.getZ(root_cs).magnitude());
-
- node["delete_on_hit"] = deleteOnHit_;
-
- return node;
-}
-
-template <typename TTracking, typename TOutput>
-inline void ObservationCubic<TTracking, TOutput>::reset() {
- energy_ = 0_GeV;
- count_ = 0;
-}
} // namespace corsika
diff --git a/corsika/detail/modules/proposal/Interaction.inl b/corsika/detail/modules/proposal/Interaction.inl
index af6d603f6874113c4430c157085557956ca19581..d8916cbe976e77b8865d7d0edaacc6bfc0135e0b 100644
--- a/corsika/detail/modules/proposal/Interaction.inl
+++ b/corsika/detail/modules/proposal/Interaction.inl
@@ -33,7 +33,7 @@ namespace corsika::proposal {
// take some minutes if you have to build the tables and cannot read the
// from disk
auto const emCut =
- calculate_kinetic_energy_threshold(code) +
+ get_kinetic_energy_threshold(code) +
get_mass(code); //! energy thresholds globally defined for individual particles
auto c = p_cross->second(media.at(comp.getHash()), emCut);
diff --git a/corsika/detail/modules/tracking/TrackingStraight.inl b/corsika/detail/modules/tracking/TrackingStraight.inl
index 545ee6c8c3eb88f1318f5262ddd7c81a0757d7af..a54155f6f3d42c6f78471180f90b07c59ad943ca 100644
--- a/corsika/detail/modules/tracking/TrackingStraight.inl
+++ b/corsika/detail/modules/tracking/TrackingStraight.inl
@@ -24,148 +24,140 @@
namespace corsika::tracking_line {
-template <typename TParticle>
-inline auto Tracking::makeStep(TParticle const &particle,
- LengthType steplength) {
- if (particle.getMomentum().getNorm() == 0_GeV) {
- return std::make_tuple(particle.getPosition(),
- particle.getMomentum() / 1_GeV);
- } // charge of the particle
- DirectionVector const dir = particle.getDirection();
- return std::make_tuple(particle.getPosition() + dir * steplength,
- dir.normalized());
-}
-
-template <typename TParticle>
-inline auto Tracking::getTrack(TParticle const &particle) {
- VelocityVector const initialVelocity =
- particle.getMomentum() / particle.getEnergy() * constants::c;
-
- auto const &initialPosition = particle.getPosition();
- CORSIKA_LOG_DEBUG("TrackingStraight pid: {}"
- " , E = {} GeV \n"
- "\tTracking pos: {} \n"
- "\tTracking p: {} GeV \n"
- "\tTracking v: {}",
- particle.getPID(), particle.getEnergy() / 1_GeV,
- initialPosition.getCoordinates(),
- particle.getMomentum().getComponents() / 1_GeV,
- initialVelocity.getComponents());
-
- // traverse the environment volume tree and find next
- // intersection
- auto [minTime, minNode] = nextIntersect(particle);
-
- return std::make_tuple(
- StraightTrajectory(Line(initialPosition, initialVelocity),
- minTime), // trajectory
- minNode); // next volume node
-}
-
-template <typename TParticle>
-inline Intersections Tracking::intersect(TParticle const &particle,
- Sphere const &sphere) {
- auto const &position = particle.getPosition();
- auto const delta = position - sphere.getCenter();
- auto const velocity =
- particle.getMomentum() / particle.getEnergy() * constants::c;
- auto const vSqNorm = velocity.getSquaredNorm();
- auto const R = sphere.getRadius();
-
- auto const vDotDelta = velocity.dot(delta);
- auto const discriminant =
- vDotDelta * vDotDelta - vSqNorm * (delta.getSquaredNorm() - R * R);
-
- if (discriminant.magnitude() > 0) {
- auto const sqDisc = sqrt(discriminant);
- auto const invDenom = 1 / vSqNorm;
-
- CORSIKA_LOG_TRACE("numericallyInside={}", sphere.contains(position));
- return Intersections((-vDotDelta - sqDisc) * invDenom,
- (-vDotDelta + sqDisc) * invDenom);
+ template <typename TParticle>
+ inline auto Tracking::makeStep(TParticle const& particle, LengthType steplength) {
+ if (particle.getMomentum().getNorm() == 0_GeV) {
+ return std::make_tuple(particle.getPosition(), particle.getMomentum() / 1_GeV);
+ } // charge of the particle
+ DirectionVector const dir = particle.getDirection();
+ return std::make_tuple(particle.getPosition() + dir * steplength, dir.normalized());
}
- return Intersections();
-}
-
-template <typename TParticle>
-inline Intersections Tracking::intersect(TParticle const &particle,
- Cubic const &cubic) {
- Point const &position = particle.getPosition();
- VelocityVector const velocity =
- particle.getMomentum() / particle.getEnergy() * constants::c;
- CoordinateSystemPtr const &cs = cubic.getCoordinateSystem();
- LengthType x0 = position.getX(cs);
- LengthType y0 = position.getY(cs);
- LengthType z0 = position.getZ(cs);
- SpeedType vx = velocity.getX(cs);
- SpeedType vy = velocity.getY(cs);
- SpeedType vz = velocity.getZ(cs);
- CORSIKA_LOG_TRACE("particle in cubic coordinate: position: ({:.3f}, {:.3f}, "
- "{:.3f}) m, veolocity: ({:.3f}, {:.3f}, {:.3f}) m/ns",
- x0 / 1_m, y0 / 1_m, z0 / 1_m, vx / (1_m / 1_ns),
- vy / (1_m / 1_ns), vz / (1_m / 1_ns));
-
- auto get_intersect_min_max = [](LengthType x0, SpeedType v0, LengthType dx) {
- auto t1 = (dx - x0) / v0;
- auto t2 = (-dx - x0) / v0;
- if (t1 > t2)
- return std::make_pair(t1, t2);
- else
- return std::make_pair(t2, t1);
- };
-
- auto [tx_max, tx_min] = get_intersect_min_max(x0, vx, cubic.getX());
- auto [ty_max, ty_min] = get_intersect_min_max(y0, vy, cubic.getY());
- auto [tz_max, tz_min] = get_intersect_min_max(z0, vz, cubic.getZ());
-
- TimeType t_exit = std::min(std::min(tx_max, ty_max), tz_max);
- TimeType t_enter = std::max(std::max(tx_min, ty_min), tz_min);
-
- CORSIKA_LOG_DEBUG("t_enter: {} ns, t_exit: {} ns", t_enter / 1_ns,
- t_exit / 1_ns);
- if ((t_exit > t_enter)) {
- if (t_enter < 0_s && t_exit > 0_s)
- CORSIKA_LOG_DEBUG("numericallyInside={}", cubic.contains(position));
- else if (t_enter < 0_s && t_exit < 0_s)
- CORSIKA_LOG_DEBUG("oppisite direction");
- return Intersections(std::move(t_enter), std::move(t_exit));
- } else
+
+ template <typename TParticle>
+ inline auto Tracking::getTrack(TParticle const& particle) {
+ VelocityVector const initialVelocity =
+ particle.getMomentum() / particle.getEnergy() * constants::c;
+
+ auto const& initialPosition = particle.getPosition();
+ CORSIKA_LOG_DEBUG(
+ "TrackingStraight pid: {}"
+ " , E = {} GeV \n"
+ "\tTracking pos: {} \n"
+ "\tTracking p: {} GeV \n"
+ "\tTracking v: {}",
+ particle.getPID(), particle.getEnergy() / 1_GeV, initialPosition.getCoordinates(),
+ particle.getMomentum().getComponents() / 1_GeV, initialVelocity.getComponents());
+
+ // traverse the environment volume tree and find next
+ // intersection
+ auto [minTime, minNode] = nextIntersect(particle);
+
+ return std::make_tuple(StraightTrajectory(Line(initialPosition, initialVelocity),
+ minTime), // trajectory
+ minNode); // next volume node
+ }
+
+ template <typename TParticle>
+ inline Intersections Tracking::intersect(TParticle const& particle,
+ Sphere const& sphere) {
+ auto const& position = particle.getPosition();
+ auto const delta = position - sphere.getCenter();
+ auto const velocity = particle.getMomentum() / particle.getEnergy() * constants::c;
+ auto const vSqNorm = velocity.getSquaredNorm();
+ auto const R = sphere.getRadius();
+
+ auto const vDotDelta = velocity.dot(delta);
+ auto const discriminant =
+ vDotDelta * vDotDelta - vSqNorm * (delta.getSquaredNorm() - R * R);
+
+ if (discriminant.magnitude() > 0) {
+ auto const sqDisc = sqrt(discriminant);
+ auto const invDenom = 1 / vSqNorm;
+
+ CORSIKA_LOG_TRACE("numericallyInside={}", sphere.contains(position));
+ return Intersections((-vDotDelta - sqDisc) * invDenom,
+ (-vDotDelta + sqDisc) * invDenom);
+ }
return Intersections();
-}
-
-template <typename TParticle, typename TBaseNodeType>
-inline Intersections Tracking::intersect(TParticle const &particle,
- TBaseNodeType const &volumeNode) {
- if (Sphere const *sphere =
- dynamic_cast<Sphere const *>(&volumeNode.getVolume());
- sphere) {
- return Tracking::intersect<TParticle>(particle, *sphere);
- } else if (Cubic const *cubic =
- dynamic_cast<Cubic const *>(&volumeNode.getVolume());
- cubic) {
- return Tracking::intersect<TParticle>(particle, *cubic);
- } else {
- throw std::runtime_error("The Volume type provided is not supported in "
- "Intersect(particle, node)");
}
-}
-template <typename TParticle>
-inline Intersections Tracking::intersect(TParticle const &particle,
- Plane const &plane) {
- auto const delta = plane.getCenter() - particle.getPosition();
- auto const velocity =
- particle.getMomentum() / particle.getEnergy() * constants::c;
- auto const n = plane.getNormal();
- auto const n_dot_v = n.dot(velocity);
+ template <typename TParticle>
+ inline Intersections Tracking::intersect(TParticle const& particle,
+ Cubic const& cubic) {
+ Point const& position = particle.getPosition();
+ VelocityVector const velocity =
+ particle.getMomentum() / particle.getEnergy() * constants::c;
+ CoordinateSystemPtr const& cs = cubic.getCoordinateSystem();
+ LengthType x0 = position.getX(cs);
+ LengthType y0 = position.getY(cs);
+ LengthType z0 = position.getZ(cs);
+ SpeedType vx = velocity.getX(cs);
+ SpeedType vy = velocity.getY(cs);
+ SpeedType vz = velocity.getZ(cs);
+ CORSIKA_LOG_TRACE(
+ "particle in cubic coordinate: position: ({:.3f}, {:.3f}, "
+ "{:.3f}) m, veolocity: ({:.3f}, {:.3f}, {:.3f}) m/ns",
+ x0 / 1_m, y0 / 1_m, z0 / 1_m, vx / (1_m / 1_ns), vy / (1_m / 1_ns),
+ vz / (1_m / 1_ns));
+
+ auto get_intersect_min_max = [](LengthType x0, SpeedType v0, LengthType dx) {
+ auto t1 = (dx - x0) / v0;
+ auto t2 = (-dx - x0) / v0;
+ if (t1 > t2)
+ return std::make_pair(t1, t2);
+ else
+ return std::make_pair(t2, t1);
+ };
+
+ auto [tx_max, tx_min] = get_intersect_min_max(x0, vx, cubic.getX());
+ auto [ty_max, ty_min] = get_intersect_min_max(y0, vy, cubic.getY());
+ auto [tz_max, tz_min] = get_intersect_min_max(z0, vz, cubic.getZ());
+
+ TimeType t_exit = std::min(std::min(tx_max, ty_max), tz_max);
+ TimeType t_enter = std::max(std::max(tx_min, ty_min), tz_min);
+
+ CORSIKA_LOG_DEBUG("t_enter: {} ns, t_exit: {} ns", t_enter / 1_ns, t_exit / 1_ns);
+ if ((t_exit > t_enter)) {
+ if (t_enter < 0_s && t_exit > 0_s)
+ CORSIKA_LOG_DEBUG("numericallyInside={}", cubic.contains(position));
+ else if (t_enter < 0_s && t_exit < 0_s)
+ CORSIKA_LOG_DEBUG("oppisite direction");
+ return Intersections(std::move(t_enter), std::move(t_exit));
+ } else
+ return Intersections();
+ }
- CORSIKA_LOG_TRACE("n_dot_v={}, delta={}, momentum={}", n_dot_v, delta,
- particle.getMomentum());
+ template <typename TParticle, typename TBaseNodeType>
+ inline Intersections Tracking::intersect(TParticle const& particle,
+ TBaseNodeType const& volumeNode) {
+ if (Sphere const* sphere = dynamic_cast<Sphere const*>(&volumeNode.getVolume());
+ sphere) {
+ return Tracking::intersect<TParticle>(particle, *sphere);
+ } else if (Cubic const* cubic = dynamic_cast<Cubic const*>(&volumeNode.getVolume());
+ cubic) {
+ return Tracking::intersect<TParticle>(particle, *cubic);
+ } else {
+ throw std::runtime_error(
+ "The Volume type provided is not supported in "
+ "Intersect(particle, node)");
+ }
+ }
- if (n_dot_v.magnitude() == 0)
- return Intersections();
- else
- return Intersections(n.dot(delta) / n_dot_v);
-}
+ template <typename TParticle>
+ inline Intersections Tracking::intersect(TParticle const& particle,
+ Plane const& plane) {
+ auto const delta = plane.getCenter() - particle.getPosition();
+ auto const velocity = particle.getMomentum() / particle.getEnergy() * constants::c;
+ auto const n = plane.getNormal();
+ auto const n_dot_v = n.dot(velocity);
+
+ CORSIKA_LOG_TRACE("n_dot_v={}, delta={}, momentum={}", n_dot_v, delta,
+ particle.getMomentum());
+
+ if (n_dot_v.magnitude() == 0)
+ return Intersections();
+ else
+ return Intersections(n.dot(delta) / n_dot_v);
+ }
} // namespace corsika::tracking_line
diff --git a/corsika/detail/modules/writers/ObservationCubicWriterParquet.inl b/corsika/detail/modules/writers/ObservationCubicWriterParquet.inl
index bd85fd0dfeffd0282c4142e2da124da9654e23a9..1a79a4b0d851ee7bb1b507e55bcd1d13d57c1b0c 100644
--- a/corsika/detail/modules/writers/ObservationCubicWriterParquet.inl
+++ b/corsika/detail/modules/writers/ObservationCubicWriterParquet.inl
@@ -50,15 +50,16 @@ namespace corsika {
inline void ObservationCubicWriterParquet::endOfLibrary() { output_.closeStreamer(); }
- inline void ObservationCubicWriterParquet::write(Code const& pid, HEPEnergyType const& energy,
- LengthType const& x, LengthType const& y, LengthType const& z,
- double nx, double ny, double nz,
- TimeType const& t) {
+ inline void ObservationCubicWriterParquet::write(
+ Code const& pid, HEPEnergyType const& energy, LengthType const& x,
+ LengthType const& y, LengthType const& z, double nx, double ny, double nz,
+ TimeType const& t) {
// write the next row - we must write `shower_` first.
*(output_.getWriter()) << shower_ << static_cast<int>(get_PDG(pid))
<< static_cast<float>(energy / 1_GeV)
- << static_cast<float>(x / 1_m) << static_cast<float>(y / 1_m) << static_cast<float>(z / 1_m)
- << static_cast<float>(nx) << static_cast<float>(ny) << static_cast<float>(nz)
+ << static_cast<float>(x / 1_m) << static_cast<float>(y / 1_m)
+ << static_cast<float>(z / 1_m) << static_cast<float>(nx)
+ << static_cast<float>(ny) << static_cast<float>(nz)
<< static_cast<float>(t / 1_ns) << parquet::EndRow;
}
diff --git a/corsika/framework/geometry/Cubic.hpp b/corsika/framework/geometry/Cubic.hpp
index 77b49a00eefc44aa710761a298e11d4e360404fb..d0aee3f1cd221032ea50ee82cab0db45336f0df6 100644
--- a/corsika/framework/geometry/Cubic.hpp
+++ b/corsika/framework/geometry/Cubic.hpp
@@ -10,16 +10,20 @@ namespace corsika {
public:
// a CoordinateSystemPtr to specify the orintation of coordinate
- Cubic(Point const& center, CoordinateSystemPtr cs, LengthType const x, LengthType const y, LengthType const z)
+ Cubic(Point const& center, CoordinateSystemPtr cs, LengthType const x,
+ LengthType const y, LengthType const z)
: center_(center)
, cs_(make_translation(cs, center.getCoordinates(cs)))
- , x_(x), y_(y), z_(z) {}
+ , x_(x)
+ , y_(y)
+ , z_(z) {}
Cubic(Point const& center, CoordinateSystemPtr cs, LengthType const side)
: center_(center)
, cs_(make_translation(cs, center.getCoordinates(cs)))
- , x_(side/2), y_(side/2), z_(side/2) {}
-
+ , x_(side / 2)
+ , y_(side / 2)
+ , z_(side / 2) {}
//! returns true if the Point p is within the sphere
bool contains(Point const& p) const override;
@@ -34,7 +38,8 @@ namespace corsika {
protected:
Point center_;
- CoordinateSystemPtr cs_; // local coordinate system with center_ in coordinate (0, 0, 0) and user defined orientation
+ CoordinateSystemPtr cs_; // local coordinate system with center_ in coordinate (0, 0,
+ // 0) and user defined orientation
LengthType x_;
LengthType y_;
LengthType z_;
diff --git a/corsika/framework/stack/SecondaryView.hpp b/corsika/framework/stack/SecondaryView.hpp
index 9ffaf150582b081cacff72a4d021403353b29dcb..a51538d583d91bd9dbbb5a962be3fad07555aeeb 100644
--- a/corsika/framework/stack/SecondaryView.hpp
+++ b/corsika/framework/stack/SecondaryView.hpp
@@ -88,9 +88,9 @@ namespace corsika {
inner_stack_value_type>
stack_value_iterator;
- typedef ConstStackIteratorInterface<
- typename std::remove_reference<TStackDataType>::type, TParticleInterface,
- MSecondaryProducer, inner_stack_value_type>
+ typedef ConstStackIteratorInterface<std::remove_reference_t<TStackDataType>,
+ TParticleInterface, MSecondaryProducer,
+ inner_stack_value_type>
const_stack_value_iterator;
/// @}
@@ -98,9 +98,8 @@ namespace corsika {
TParticleInterface, MSecondaryProducer, view_type>
stack_view_iterator;
- typedef ConstStackIteratorInterface<
- typename std::remove_reference<TStackDataType>::type, TParticleInterface,
- MSecondaryProducer, view_type>
+ typedef ConstStackIteratorInterface<std::remove_reference_t<TStackDataType>,
+ TParticleInterface, MSecondaryProducer, view_type>
const_stack_view_iterator;
/**
diff --git a/corsika/modules/ObservationCubic.hpp b/corsika/modules/ObservationCubic.hpp
index f38d23d6775246d8e103d944c0e78c0fab5d8e58..f7179bf117f0651126634212340e3ab29d80a55d 100644
--- a/corsika/modules/ObservationCubic.hpp
+++ b/corsika/modules/ObservationCubic.hpp
@@ -5,7 +5,6 @@
#include <corsika/modules/writers/ObservationCubicWriterParquet.hpp>
#include <corsika/framework/process/ContinuousProcess.hpp>
-
namespace corsika {
/**
@@ -30,8 +29,8 @@ namespace corsika {
public TOutputWriter {
public:
- ObservationCubic(Point const& center, CoordinateSystemPtr cs,
- LengthType const x, LengthType const y, LengthType const z, bool = true);
+ ObservationCubic(Point const& center, CoordinateSystemPtr cs, LengthType const x,
+ LengthType const y, LengthType const z, bool = true);
template <typename TParticle, typename TTrajectory>
ProcessReturn doContinuous(TParticle& vParticle, TTrajectory& vTrajectory,
@@ -49,11 +48,8 @@ namespace corsika {
bool const deleteOnHit_;
HEPEnergyType energy_;
unsigned int count_;
-
};
//! @}
} // namespace corsika
-
-
#include <corsika/detail/modules/ObservationCubic.inl>
\ No newline at end of file
diff --git a/corsika/modules/writers/ObservationCubicWriterParquet.hpp b/corsika/modules/writers/ObservationCubicWriterParquet.hpp
index 449cb380380aa3cd762ee7d0e9bb20846a4bdf15..97e518c8d3da9e9272dce572cf46cbfdc4429d55 100644
--- a/corsika/modules/writers/ObservationCubicWriterParquet.hpp
+++ b/corsika/modules/writers/ObservationCubicWriterParquet.hpp
@@ -49,9 +49,8 @@ namespace corsika {
/**
* Write a particle to the file.
*/
- void write(Code const& pid, HEPEnergyType const& energy,
- LengthType const& x, LengthType const& y, LengthType const& z,
- double nx, double ny, double nz,
+ void write(Code const& pid, HEPEnergyType const& energy, LengthType const& x,
+ LengthType const& y, LengthType const& z, double nx, double ny, double nz,
TimeType const& t);
}; // class ObservationCubicWriterParquet