From cc13c82c1ae0d79fcb944f1a3181b7fb07f97c87 Mon Sep 17 00:00:00 2001
From: Andre Schmidt <schmidt-a@iklx288.ikp.kit.edu>
Date: Tue, 13 Oct 2020 17:50:43 +0200
Subject: [PATCH] min distance is no longer steplength
---
Framework/Cascade/Cascade.h | 17 +++++++++++------
Processes/ObservationPlane/ObservationPlane.cc | 2 +-
2 files changed, 12 insertions(+), 7 deletions(-)
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index 1b41d0a68..9896259cb 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -408,7 +408,7 @@ namespace corsika::cascade {
// determine steplength for the magnetic field
// because Steplength should not be min_distance
- /*
+
int chargeNumber;
if (corsika::particles::IsNucleus(vParticle.GetPID())) {
chargeNumber = vParticle.GetNuclearZ();
@@ -419,16 +419,21 @@ namespace corsika::cascade {
auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(vParticle.GetPosition());
geometry::Vector<dimensionless_d> const directionBefore = vParticle.GetMomentum().normalized();
auto c = directionBefore.cross(magneticfield) * chargeNumber * corsika::units::constants::c * 1_eV /
- (vParticle.GetMomentum().norm() * 1_V)
+ (vParticle.GetMomentum().norm() * 1_V);
LengthType Steplength = min_distance;
if (chargeNumber != 0) {
Steplength = sqrt(2 / c.squaredNorm() * (sqrt(c.squaredNorm() * min_distance * min_distance + 1) -1));
std::cout << "Steplength " << Steplength << std::endl;
- }
- */
+ }
+ if (Steplength == 0_m) {
+ Steplength = min_distance;
+ }
+
// This formula hasnt been tested
- auto [position, direction, L2] = fTracking.MagneticStep(vParticle, min_distance);
+ //auto [position, direction, L2] = fTracking.MagneticStep(vParticle, min_distance);
+ auto [position, direction, L2] = fTracking.MagneticStep(vParticle, Steplength);
+
//histL2(L2);
histLlog2(L2);
int pdg = static_cast<int>(particles::GetPDG(vParticle.GetPID()));
@@ -481,7 +486,7 @@ namespace corsika::cascade {
- int chargeNumber = 0;
+ //int chargeNumber = 0;
if (corsika::particles::IsNucleus(vParticle.GetPID())) {
chargeNumber = vParticle.GetNuclearZ();
} else {
diff --git a/Processes/ObservationPlane/ObservationPlane.cc b/Processes/ObservationPlane/ObservationPlane.cc
index 036bbcacf..ff70be5fe 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());
auto direction = trajectory.GetV0().normalized();
- if (chargeNumber != 0 && plane_.GetNormal().dot(trajectory.GetV0().cross(magneticfield)) * 1_s / 1_m / 1_T > 1e-6) {
+ if (chargeNumber != 0 && abs(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::c * 1_eV / (vParticle.GetMomentum().norm() * 1_V);
--
GitLab