From f2f807faf40ff8faefc9443b9a081a75b44f18b1 Mon Sep 17 00:00:00 2001
From: Andre Schmidt <schmidt-a@ikp-admins@lists.kit.edu>
Date: Mon, 21 Sep 2020 13:44:23 +0200
Subject: [PATCH] momentum
---
.../ObservationPlane/ObservationPlane.cc | 26 ++++++++++---------
Processes/TrackingLine/TrackingLine.h | 4 +--
2 files changed, 16 insertions(+), 14 deletions(-)
diff --git a/Processes/ObservationPlane/ObservationPlane.cc b/Processes/ObservationPlane/ObservationPlane.cc
index 27bed61ac..036bbcacf 100644
--- a/Processes/ObservationPlane/ObservationPlane.cc
+++ b/Processes/ObservationPlane/ObservationPlane.cc
@@ -69,25 +69,26 @@ LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const& vPa
}
auto const* currentLogicalVolumeNode = vParticle.GetNode();
auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(vParticle.GetPosition());
- geometry::Vector<SpeedType::dimension_type> const velocity = trajectory.GetV0();
+ auto direction = trajectory.GetV0().normalized();
- if (chargeNumber != 0 && plane_.GetNormal().dot(velocity.cross(magneticfield)) * 1_s / 1_m / 1_T > 0.0001) {
+ if (chargeNumber != 0 && plane_.GetNormal().dot(trajectory.GetV0().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 * corsika::units::constants::cSquared * 1_eV /
- (velocity.GetSquaredNorm() * vParticle.GetEnergy() * 1_V);
+ auto k = chargeNumber * corsika::units::constants::c * 1_eV / (vParticle.GetMomentum().norm() * 1_V);
+
LengthType MaxStepLength1 =
- ( sqrt(velocity.dot(plane_.GetNormal()) * velocity.dot(plane_.GetNormal()) / velocity.GetSquaredNorm() -
+ ( sqrt(direction.dot(plane_.GetNormal()) * direction.dot(plane_.GetNormal()) -
(plane_.GetNormal().dot(trajectory.GetR0() - plane_.GetCenter()) *
- plane_.GetNormal().dot(velocity.cross(magneticfield)) * 2 * k)) -
- velocity.dot(plane_.GetNormal()) / velocity.GetNorm() ) /
- (plane_.GetNormal().dot(velocity.cross(magneticfield)) * k);
+ plane_.GetNormal().dot(direction.cross(magneticfield)) * 2 * k)) -
+ direction.dot(plane_.GetNormal()) / direction.GetNorm() ) /
+ (plane_.GetNormal().dot(direction.cross(magneticfield)) * k);
+
LengthType MaxStepLength2 =
- ( - sqrt(velocity.dot(plane_.GetNormal()) * velocity.dot(plane_.GetNormal()) / velocity.GetSquaredNorm() -
+ ( - sqrt(direction.dot(plane_.GetNormal()) * direction.dot(plane_.GetNormal()) -
(plane_.GetNormal().dot(trajectory.GetR0() - plane_.GetCenter()) *
- plane_.GetNormal().dot(velocity.cross(magneticfield)) * 2 * k)) -
- velocity.dot(plane_.GetNormal()) / velocity.GetNorm() ) /
- (plane_.GetNormal().dot(velocity.cross(magneticfield)) * k);
+ 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;
@@ -108,6 +109,7 @@ LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const& vPa
}
auto const pointOfIntersection = trajectory.GetPosition(timeOfIntersection);
+ std::cout << " obs plane non b-field " << std::endl;
return (trajectory.GetR0() - pointOfIntersection).norm() * 1.0001;
}
diff --git a/Processes/TrackingLine/TrackingLine.h b/Processes/TrackingLine/TrackingLine.h
index 10f848089..b1906de69 100644
--- a/Processes/TrackingLine/TrackingLine.h
+++ b/Processes/TrackingLine/TrackingLine.h
@@ -293,7 +293,7 @@ namespace corsika::process {
// creating Line with magnetic field
if (chargeNumber != 0) {
- auto k = chargeNumber * corsika::units::constants::cSquared * 1_eV / (velocity.norm() * p.GetEnergy() * 1_V);
+ auto k = chargeNumber * corsika::units::constants::c * 1_eV / (p.GetMomentum().norm() * 1_V);
geometry::Vector<dimensionless_d> const directionBefore = velocity.normalized();
// determine steplength to next volume
double a = ((directionBefore.cross(magneticfield)).dot(currentPosition - sphere.GetCenter()) * k + 1) * 4 /
@@ -352,7 +352,7 @@ namespace corsika::process {
// creating Line with magnetic field
if (chargeNumber != 0) {
- auto k = chargeNumber * corsika::units::constants::cSquared * 1_eV / (velocity.norm() * p.GetEnergy() * 1_V);
+ auto k = chargeNumber * corsika::units::constants::c * 1_eV / (p.GetMomentum().norm() * 1_V);
geometry::Vector<dimensionless_d> const directionBefore = velocity.normalized();
// determine steplength to next volume
double a = ((directionBefore.cross(magneticfield)).dot(currentPosition - sphere.GetCenter()) * k + 1) * 4 /
--
GitLab