From bb1a2819a28a588f7c295f0d6fa52b3e8db31a12 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Fri, 22 Jan 2021 11:35:18 +0000
Subject: [PATCH] ObservationPlane
---
corsika/detail/modules/ObservationPlane.inl | 69 ++++++++++++++++++---
1 file changed, 60 insertions(+), 9 deletions(-)
diff --git a/corsika/detail/modules/ObservationPlane.inl b/corsika/detail/modules/ObservationPlane.inl
index 74c29b336..5c521943c 100644
--- a/corsika/detail/modules/ObservationPlane.inl
+++ b/corsika/detail/modules/ObservationPlane.inl
@@ -59,15 +59,66 @@ namespace corsika {
corsika::setup::Stack::particle_type const& particle,
corsika::setup::Trajectory const& trajectory) {
- 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());
+ 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;
+ }
+ }
TimeType const timeOfIntersection = intersection.getEntry();
--
GitLab