adaptions

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);

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

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)