From 2c3c6d932b40cdbb9ce32ab3b8a1319ec4da46a2 Mon Sep 17 00:00:00 2001
From: Andre Schmidt <schmidt-a@ikp-admins@lists.kit.edu>
Date: Tue, 18 Aug 2020 15:31:18 +0200
Subject: [PATCH] solved transition error
---
Framework/Cascade/Cascade.h | 23 ++++++++-----------
.../ObservationPlane/ObservationPlane.cc | 5 ++--
Processes/TrackingLine/TrackingLine.h | 13 ++++++++---
3 files changed, 23 insertions(+), 18 deletions(-)
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index 9ba4c4bc8..9d0a296d5 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -229,18 +229,15 @@ namespace corsika::cascade {
// take minimum of geometry, interaction, decay for next step
auto min_distance = std::min(
{distance_interact, distance_decay, distance_max, geomMaxLength});
- if (min_distance == geomMaxLength) {
- min_distance = 1.001 * geomMaxLength;
- }
C8LOG_DEBUG("transport particle by : {} m", min_distance / 1_m);
// determine displacement by the magnetic field
- auto const* currentLogicalVolumeNode = vParticle.GetNode();
+ auto const* currentLogicalVolumeNode = vParticle.GetNode();
auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(vParticle.GetPosition());
geometry::Vector<SpeedType::dimension_type> velocity = vParticle.GetMomentum() / vParticle.GetEnergy() *
corsika::units::constants::c;
- geometry::Vector<dimensionless_d> const directionBefore = velocity.normalized();
+ geometry::Vector<dimensionless_d> const directionBefore = velocity.normalized();
int chargeNumber;
if (corsika::particles::IsNucleus(vParticle.GetPID())) {
chargeNumber = vParticle.GetNuclearZ();
@@ -250,19 +247,19 @@ namespace corsika::cascade {
auto k = chargeNumber * corsika::units::constants::cSquared * 1_eV /
(velocity.GetNorm() * vParticle.GetEnergy() * 1_V);
- // First Movement
- // assuming magnetic field does not change during movement
- auto position = vParticle.GetPosition() + directionBefore * min_distance / 2;
- // Change of direction by magnetic field
- geometry::Vector<dimensionless_d> const directionAfter = directionBefore + directionBefore.cross(magneticfield) *
+ // First Movement
+ // assuming magnetic field does not change during movement
+ auto position = vParticle.GetPosition() + directionBefore * min_distance / 2;
+ // Change of direction by magnetic field
+ geometry::Vector<dimensionless_d> const directionAfter = directionBefore + directionBefore.cross(magneticfield) *
min_distance * k;
- // Second Movement
- position = position + directionAfter * min_distance / 2;
+ // Second Movement
+ position = position + directionAfter * min_distance / 2;
auto distance = position - vParticle.GetPosition();
// distance.norm() != min_distance if q != 0
// small error can be neglected
if (distance.norm() != 0_m) {
- velocity = distance.normalized() * velocity.norm();
+ velocity = distance.normalized() * velocity.norm();
} // no velocity update for very small steps
// here the particle is actually moved along the trajectory to new position:
diff --git a/Processes/ObservationPlane/ObservationPlane.cc b/Processes/ObservationPlane/ObservationPlane.cc
index d12a55f5a..b135f9c7a 100644
--- a/Processes/ObservationPlane/ObservationPlane.cc
+++ b/Processes/ObservationPlane/ObservationPlane.cc
@@ -71,7 +71,7 @@ LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const& vPa
auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(vParticle.GetPosition());
geometry::Vector<SpeedType::dimension_type> const velocity = trajectory.GetV0();
- if (chargeNumber != 0 && plane_.GetNormal().dot(velocity.cross(magneticfield)) * 1_s / 1_m / 1_T != 0) {
+ if (chargeNumber != 0 && plane_.GetNormal().dot(velocity.cross(magneticfield)) * 1_s / 1_m / 1_T > 0.0001) {
auto const* currentLogicalVolumeNode = vParticle.GetNode();
auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(vParticle.GetPosition());
auto k = chargeNumber * corsika::units::constants::cSquared * 1_eV /
@@ -88,12 +88,13 @@ LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const& vPa
plane_.GetNormal().dot(velocity.cross(magneticfield)) * 2 * k)) -
velocity.dot(plane_.GetNormal()) / velocity.GetNorm() ) /
(plane_.GetNormal().dot(velocity.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) {
+ std::cout << " distance to obs plane : " << MaxStepLength2 << std::endl;
return MaxStepLength2 * 1.0001;
} else if (MaxStepLength2 <= 0_m || MaxStepLength1 < MaxStepLength2) {
+ std::cout << " distance to obs plane : " << MaxStepLength1 << std::endl;
return MaxStepLength1 * 1.0001;
}
}
diff --git a/Processes/TrackingLine/TrackingLine.h b/Processes/TrackingLine/TrackingLine.h
index eff54e7e2..a888cffaa 100644
--- a/Processes/TrackingLine/TrackingLine.h
+++ b/Processes/TrackingLine/TrackingLine.h
@@ -111,12 +111,13 @@ namespace corsika::process {
for (int i = 0; i < 4; i++) {
if (solutions[i].imag() == 0 && solutions[i].real() > 0) {
tmp.push_back(solutions[i].real());
+ std::cout << "Solutions for next Volume: " << solutions[i].real() << std::endl;
}
}
LengthType Steplength;
if (tmp.size() > 0) {
Steplength = 1_m * *std::min_element(tmp.begin(),tmp.end());
- std::cout << "s = " << Steplength << std::endl;
+ std::cout << "Steplength to next volume = " << Steplength << std::endl;
} else {
std::cout << "no intersection (1)!" << std::endl;
// what to do when this happens? (very unlikely)
@@ -172,12 +173,18 @@ namespace corsika::process {
for (int i = 0; i < 4; i++) {
if (solutions[i].imag() == 0 && solutions[i].real() > 0) {
tmp.push_back(solutions[i].real());
+ std::cout << "Solutions for current Volume: " << solutions[i].real() << std::endl;
}
}
LengthType Steplength;
if (tmp.size() > 0) {
- Steplength = 1_m * *std::min_element(tmp.begin(),tmp.end());
- std::cout << "s = " << Steplength << std::endl;
+ Steplength = 1_m * *std::min_element(tmp.begin(),tmp.end());
+ if (numericallyInside == false) {
+ int p = std::min_element(tmp.begin(),tmp.end()) - tmp.begin();
+ tmp.erase(tmp.begin() + p);
+ Steplength = 1_m * *std::min_element(tmp.begin(),tmp.end());
+ }
+ std::cout << "steplength out of current volume = " << Steplength << std::endl;
} else {
std::cout << "no intersection (2)!" << std::endl;
// what to do when this happens? (very unlikely)
--
GitLab