From 951cdae2c16bcec492ea56b1dafbaa4f9cf792a1 Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Thu, 28 Jan 2021 22:30:15 +0100
Subject: [PATCH] adaptions
---
corsika/detail/modules/ObservationPlane.inl | 72 ++-----------------
corsika/detail/modules/OnShellCheck.inl | 10 +--
.../detail/modules/pythia8/Interaction.inl | 4 +-
3 files changed, 13 insertions(+), 73 deletions(-)
diff --git a/corsika/detail/modules/ObservationPlane.inl b/corsika/detail/modules/ObservationPlane.inl
index 5c521943c..13f605efc 100644
--- a/corsika/detail/modules/ObservationPlane.inl
+++ b/corsika/detail/modules/ObservationPlane.inl
@@ -59,82 +59,24 @@ namespace corsika {
corsika::setup::Stack::particle_type const& particle,
corsika::setup::Trajectory const& trajectory) {
- int chargeNumber;
- if (is_nucleus(vParticle.getPID())) {
- chargeNumber = vParticle.getNuclearZ();
- } else {
- chargeNumber = get_charge_number(vParticle.getPID());
- }
- auto const* currentLogicalVolumeNode = vParticle.getNode();
- auto magneticfield = currentLogicalVolumeNode->getModelProperties().getMagneticField(
- vParticle.getPosition());
- auto direction = trajectory.getVelocity(0).normalized();
-
- if (chargeNumber != 0 &&
- abs(plane_.getNormal().dot(
- trajectory.getLine().getVelocity().cross(magneticfield))) *
- 1_s / 1_m / 1_T >
- 1e-6) {
- auto const* currentLogicalVolumeNode = vParticle.getNode();
- auto magneticfield =
- currentLogicalVolumeNode->getModelProperties().getMagneticField(
- vParticle.getPosition());
- auto k =
- chargeNumber * constants::c * 1_eV / (vParticle.getMomentum().getNorm() * 1_V);
-
- if (direction.dot(plane_.getNormal()) * direction.dot(plane_.getNormal()) -
- (plane_.getNormal().dot(trajectory.getPosition(0) - plane_.getCenter()) *
- plane_.getNormal().dot(direction.cross(magneticfield)) * 2 * k) <
- 0) {
- return std::numeric_limits<double>::infinity() * 1_m;
- }
-
- LengthType MaxStepLength1 =
- (sqrt(direction.dot(plane_.getNormal()) * direction.dot(plane_.getNormal()) -
- (plane_.getNormal().dot(trajectory.getPosition(0) - plane_.getCenter()) *
- plane_.getNormal().dot(direction.cross(magneticfield)) * 2 * k)) -
- direction.dot(plane_.getNormal()) / direction.getNorm()) /
- (plane_.getNormal().dot(direction.cross(magneticfield)) * k);
-
- LengthType MaxStepLength2 =
- (-sqrt(direction.dot(plane_.getNormal()) * direction.dot(plane_.getNormal()) -
- (plane_.getNormal().dot(trajectory.getPosition(0) - plane_.getCenter()) *
- plane_.getNormal().dot(direction.cross(magneticfield)) * 2 * k)) -
- direction.dot(plane_.getNormal()) / direction.getNorm()) /
- (plane_.getNormal().dot(direction.cross(magneticfield)) * k);
-
- if (MaxStepLength1 <= 0_m && MaxStepLength2 <= 0_m) {
- return std::numeric_limits<double>::infinity() * 1_m;
- } else if (MaxStepLength1 <= 0_m || MaxStepLength2 < MaxStepLength1) {
- CORSIKA_LOG_DEBUG("steplength to obs plane 2: {} ", MaxStepLength2);
- return MaxStepLength2 *
- (direction + direction.cross(magneticfield) * MaxStepLength2 * k / 2)
- .getNorm() *
- 1.001;
- } else if (MaxStepLength2 <= 0_m || MaxStepLength1 < MaxStepLength2) {
- CORSIKA_LOG_DEBUG("steplength to obs plane 1: {}", MaxStepLength1);
- return MaxStepLength1 *
- (direction + direction.cross(magneticfield) * MaxStepLength2 * k / 2)
- .getNorm() *
- 1.001;
- }
- }
-
+ auto const& volumeNode = particle.getNode();
+ typedef typename std::remove_const_t<
+ std::remove_reference_t<decltype(volumeNode->getModelProperties())>>
+ medium_type;
+ Intersections const intersection =
+ setup::Tracking::intersect<corsika::setup::Stack::particle_type, medium_type>(
+ particle, plane_, volumeNode->getModelProperties());
TimeType const timeOfIntersection = intersection.getEntry();
-
CORSIKA_LOG_TRACE("particle={}, pos={}, dir={}, plane={}, timeOfIntersection={}",
particle.asString(), particle.getPosition(),
particle.getDirection(), plane_.asString(), 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();
-
auto const pointOfIntersection = trajectory.getPosition(fractionOfIntersection);
auto dist = (trajectory.getPosition(0) - pointOfIntersection).getNorm();
CORSIKA_LOG_TRACE("ObservationPlane: getMaxStepLength l={} m", dist / 1_m);
diff --git a/corsika/detail/modules/OnShellCheck.inl b/corsika/detail/modules/OnShellCheck.inl
index bd77def88..2702de6e4 100644
--- a/corsika/detail/modules/OnShellCheck.inl
+++ b/corsika/detail/modules/OnShellCheck.inl
@@ -73,14 +73,14 @@ namespace corsika {
if (throw_error_)
throw std::runtime_error(
"OnShellCheck: error! shifted energy by large amount!");
- }
}
+ }
- // reset energy
- p.setEnergy(e_shifted);
- } else
- CORSIKA_LOGGER_DEBUG(logger_, "particle mass for {} OK", pid);
+ // reset energy
+ p.setEnergy(e_shifted);
}
+ else CORSIKA_LOGGER_DEBUG(logger_, "particle mass for {} OK", pid);
}
+} // namespace corsika
} // namespace corsika
diff --git a/corsika/detail/modules/pythia8/Interaction.inl b/corsika/detail/modules/pythia8/Interaction.inl
index 8167639e5..eb30b23a6 100644
--- a/corsika/detail/modules/pythia8/Interaction.inl
+++ b/corsika/detail/modules/pythia8/Interaction.inl
@@ -21,9 +21,7 @@
namespace corsika::pythia8 {
- Interaction::~Interaction() {
- CORSIKA_LOG_INFO("Pythia::Interaction n= {}", count_);
- }
+ Interaction::~Interaction() { CORSIKA_LOG_INFO("Pythia::Interaction n= {}", count_); }
Interaction::Interaction(bool const print_listing)
: Pythia8::Pythia(CORSIKA_Pythia8_XML_DIR)
--
GitLab