From eb14f7fd130339cfff7f14b9f44800823c5b1311 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Andr=C3=A9=20Schmidt?= <upwli@student.kit.edu>
Date: Thu, 9 Jul 2020 11:49:18 +0200
Subject: [PATCH] Replace Cascade.h
---
Framework/Cascade/Cascade.h | 53 ++++++++++++++++++++-----------------
1 file changed, 29 insertions(+), 24 deletions(-)
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index 500443b12..8779ed513 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -220,31 +220,36 @@ namespace corsika::cascade {
} else {
chargeNumber = corsika::particles::GetChargeNumber(vParticle.GetPID());
}
- auto magneticfield = corsika::geometry::Vector
- (fEnvironment.GetCoordinateSystem(), 0_uT, 50_uT, 0_uT);
geometry::Vector<SpeedType::dimension_type> const velocity =
- vParticle.GetMomentum() / vParticle.GetEnergy() * corsika::units::constants::c;
- geometry::Vector<SpeedType::dimension_type> const velocityVerticalMag =
- velocity - magneticfield * velocity.dot(magneticfield) / (magneticfield.GetSquaredNorm());
- geometry::Vector<dimensionless_d> const directionBefore = velocity / velocity.GetNorm();
- LengthType const gyroradius = vParticle.GetEnergy() * velocityVerticalMag.GetNorm() * 1_V /
- (corsika::units::constants::cSquared * abs(chargeNumber) *
- magneticfield.GetNorm() * 1_eV);
- LengthType const Steplength = 0.01 * gyroradius;
- // First Movement
- auto position = vParticle.GetPosition() + directionBefore * Steplength / 2;
- // Change of direction by magnetic field at position
- magneticfield = corsika::geometry::Vector(fEnvironment.GetCoordinateSystem(), 0_uT, 50_uT, 0_uT);
- geometry::Vector<dimensionless_d> const directionAfter = directionBefore +
- velocity.cross(magneticfield) * chargeNumber *
- Steplength * corsika::units::constants::cSquared * 1_eV/
- (vParticle.GetEnergy() * velocity.GetSquaredNorm() * 1_V);
- // Second Movement
- position = position + directionAfter * Steplength / 2;
- LengthType const magMaxLength = (position - vParticle.GetPosition()).GetNorm();
- geometry::Vector<dimensionless_d> const direction = (position - vParticle.GetPosition()) /
- magMaxLength;
- vParticle.SetMomentum( direction * vParticle.GetMomentum().GetNorm());
+ vParticle.GetMomentum() / vParticle.GetEnergy() * corsika::units::constants::c;
+ geometry::Vector<dimensionless_d> const directionBefore = velocity / velocity.GetNorm();
+ auto magMaxLength = 1_m/0;
+ auto directionAfter = directionBefore;
+ if(chargeNumber != 0) {
+ auto magneticfield = corsika::geometry::Vector
+ (fEnvironment.GetCoordinateSystem(), 0_uT, 50_uT, 0_uT);
+ geometry::Vector<SpeedType::dimension_type> const velocityVerticalMag =
+ velocity - magneticfield * velocity.dot(magneticfield) / (magneticfield.GetSquaredNorm());
+ LengthType const gyroradius = vParticle.GetEnergy() * velocityVerticalMag.GetNorm() * 1_V /
+ (corsika::units::constants::cSquared * abs(chargeNumber) *
+ magneticfield.GetNorm() * 1_eV);
+ //steplength depending on how exact it should
+ LengthType const Steplength = 0.01 * gyroradius;
+ // First Movement
+ auto position = vParticle.GetPosition() + directionBefore * Steplength / 2;
+ // Change of direction by magnetic field at position
+ magneticfield = corsika::geometry::Vector(fEnvironment.GetCoordinateSystem(), 0_uT, 50_uT, 0_uT);
+ directionAfter = directionBefore + velocity.cross(magneticfield) * chargeNumber *
+ Steplength * corsika::units::constants::cSquared * 1_eV /
+ (vParticle.GetEnergy() * velocity.GetSquaredNorm() * 1_V);
+ // Second Movement
+ position = position + directionAfter * Steplength / 2;
+ magMaxLength = (position - vParticle.GetPosition()).GetNorm();
+ geometry::Vector<dimensionless_d> const direction = (position - vParticle.GetPosition()) /
+ magMaxLength;
+ vParticle.SetMomentum( direction * vParticle.GetMomentum().GetNorm());
+ std::cout << "New direction because of magnetic field: " << direction.GetComponents() << std::endl;
+ }
// determine geometric tracking
auto [step, geomMaxLength, nextVol] = fTracking.GetTrack(vParticle);
--
GitLab