From d3c475e5f9abc2cc4d3374dc87d6dc1fea49a9a9 Mon Sep 17 00:00:00 2001
From: Andre Schmidt <schmidt-a@iklx288.ikp.kit.edu>
Date: Sat, 22 Aug 2020 08:55:50 +0200
Subject: [PATCH] added function for magnetic step
---
Framework/Cascade/Cascade.h | 33 +-------
Processes/TrackingLine/TrackingLine.h | 117 ++++++++++++++------------
2 files changed, 66 insertions(+), 84 deletions(-)
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index 9d0a296d5..f4b8e92d2 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -214,7 +214,7 @@ namespace corsika::cascade {
vParticle.GetEnergy() * units::constants::c;
// determine geometric tracking
- auto [stepWithoutB, stepWithB, geomMaxLength, nextVol] = fTracking.GetTrack(vParticle);
+ auto [lineWithoutB, stepWithoutB, stepWithB, geomMaxLength, nextVol] = fTracking.GetTrack(vParticle);
[[maybe_unused]] auto const& dummy_nextVol = nextVol;
// convert next_step from grammage to length
@@ -232,41 +232,16 @@ namespace corsika::cascade {
C8LOG_DEBUG("transport particle by : {} m", min_distance / 1_m);
- // determine displacement by the magnetic field
- 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();
- int chargeNumber;
- if (corsika::particles::IsNucleus(vParticle.GetPID())) {
- chargeNumber = vParticle.GetNuclearZ();
- } else {
- chargeNumber = corsika::particles::GetChargeNumber(vParticle.GetPID());
- }
- 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) *
- min_distance * k;
- // Second Movement
- position = position + directionAfter * min_distance / 2;
+ // determine displacement by the magnetic field
+ auto [line, position, directionAfter] = fTracking.MagneticStep(vParticle, lineWithoutB, min_distance);
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();
- } // no velocity update for very small steps
// here the particle is actually moved along the trajectory to new position:
// std::visit(setup::ParticleUpdate<Particle>{vParticle}, step);
vParticle.SetMomentum(directionAfter.normalized() * vParticle.GetMomentum().norm());
- geometry::Line line(vParticle.GetPosition(), velocity);
- geometry::Trajectory<geometry::Line> stepNew(line, distance.norm() / velocity.norm());
+ geometry::Trajectory<geometry::Line> stepNew(line, distance.norm() / line.GetV0().norm());
vParticle.SetPosition(position);
vParticle.SetTime(vParticle.GetTime() + distance.norm() / units::constants::c);
std::cout << "New Position: " << vParticle.GetPosition().GetCoordinates() << std::endl;
diff --git a/Processes/TrackingLine/TrackingLine.h b/Processes/TrackingLine/TrackingLine.h
index a888cffaa..d68a96d65 100644
--- a/Processes/TrackingLine/TrackingLine.h
+++ b/Processes/TrackingLine/TrackingLine.h
@@ -62,6 +62,8 @@ namespace corsika::process {
<< " GeV " << std::endl;
std::cout << "TrackingLine v: " << velocity.GetComponents() << std::endl;
+ geometry::Line line(currentPosition, velocity);
+
auto const* currentLogicalVolumeNode = p.GetNode();
//~ auto const* currentNumericalVolumeNode =
//~ fEnvironment.GetUniverse()->GetContainingNode(currentPosition);
@@ -86,8 +88,7 @@ namespace corsika::process {
std::cout << " Magnetic Field: " << magneticfield.GetComponents() / 1_uT << " uT " << std::endl;
auto k = chargeNumber * corsika::units::constants::cSquared * 1_eV / (velocity.GetNorm() * p.GetEnergy() * 1_V);
geometry::Vector<dimensionless_d> const directionBefore = velocity.normalized();
- geometry::Vector<SpeedType::dimension_type> velocity1 = velocity;
- geometry::Vector<SpeedType::dimension_type> velocity2 = velocity;
+ LengthType Steplength = 10_m; // length irrelevant if q=0 and else it gets changed again
// for entering from outside
auto addIfIntersects = [&](auto const& vtn) {
@@ -114,7 +115,6 @@ namespace corsika::process {
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 << "Steplength to next volume = " << Steplength << std::endl;
@@ -123,23 +123,15 @@ namespace corsika::process {
// what to do when this happens? (very unlikely)
}
delete [] solutions;
-
- // First Movement
- // assuming magnetic field does not change during movement
- auto position = currentPosition + directionBefore * Steplength / 2;
- // Change of direction by magnetic field
- geometry::Vector<dimensionless_d> const directionAfter = directionBefore + directionBefore.cross(magneticfield) *
- Steplength * k;
- // Second Movement
- position = position + directionAfter * Steplength / 2;
- geometry::Vector<dimensionless_d> const direction = (position - currentPosition) /
- (position - currentPosition).GetNorm();
- velocity1 = direction * velocity.GetNorm();
- } // instead of changing the line with magnetic field, the TimeOfIntersection() could be changed
+
+ }
+
+ auto [line1, position, direction] = MagneticStep(p, line, Steplength);
+ // new particle position and direction not needed in this case
+ // instead of changing the line with magnetic field, the TimeOfIntersection() could be changed
// using line has some errors for huge steps
- geometry::Line line(currentPosition, velocity1);
- if (auto opt = TimeOfIntersection(line, sphere); opt.has_value()) {
+ if (auto opt = TimeOfIntersection(line1, sphere); opt.has_value()) {
auto const [t1, t2] = *opt;
C8LOG_DEBUG("intersection times: {} s; {} s", t1 / 1_s, t2 / 1_s);
if (t1.magnitude() > 0)
@@ -178,34 +170,26 @@ namespace corsika::process {
}
LengthType Steplength;
if (tmp.size() > 0) {
- 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;
+ 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)
}
delete [] solutions;
+ }
- // First Movement
- // assuming magnetic field does not change during movement
- auto position = currentPosition + directionBefore * Steplength / 2;
- // Change of direction by magnetic field
- geometry::Vector<dimensionless_d> const directionAfter = directionBefore + directionBefore.cross(magneticfield) *
- Steplength * k;
- // Second Movement
- position = position + directionAfter * Steplength / 2;
- geometry::Vector<dimensionless_d> const direction = (position - currentPosition) /
- (position - currentPosition).GetNorm();
- velocity2 = direction * velocity.GetNorm();
- } // instead of changing the line with magnetic field, the TimeOfIntersection() could be changed
- geometry::Line line(currentPosition, velocity2);
+ auto [line2, position, direction] = MagneticStep(p, line, Steplength);
+ // new particle position and direction not needed in this case
+ // instead of changing the line with magnetic field, the TimeOfIntersection() could be changed
+ // using line has some errors for huge steps
- [[maybe_unused]] auto const [t1, t2] = *TimeOfIntersection(line, sphere);
+ [[maybe_unused]] auto const [t1, t2] = *TimeOfIntersection(line2, sphere);
[[maybe_unused]] auto dummy_t1 = t1;
intersections.emplace_back(t2, currentLogicalVolumeNode->GetParent());
}
@@ -228,26 +212,49 @@ namespace corsika::process {
<< min
// << " " << minIter->second->GetModelProperties().GetName()
<< std::endl;
+
+ auto [lineWithB, position, direction] = MagneticStep(p, line, velocity.norm() * min);
+ // new particle position and direction not needed in this case
+
+ return std::make_tuple(line, geometry::Trajectory<geometry::Line>(line, min),
+ geometry::Trajectory<geometry::Line>(lineWithB, min),
+ velocity.norm() * min, minIter->second);
+ }
+
+ template <typename Particle> // was Stack previously, and argument was
+ // Stack::StackIterator
+
+ // determine direction of the particle after adding magnetic field
+ auto MagneticStep(Particle const& p, corsika::geometry::Line line, corsika::units::si::LengthType Steplength) {
+ using namespace corsika::units::si;
+
+ //charge of the particle
+ int chargeNumber;
+ if (corsika::particles::IsNucleus(p.GetPID())) {
+ chargeNumber = p.GetNuclearZ();
+ } else {
+ chargeNumber = corsika::particles::GetChargeNumber(p.GetPID());
+ }
+ auto const* currentLogicalVolumeNode = p.GetNode();
+ auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(p.GetPosition());
+ auto k = chargeNumber * corsika::units::constants::cSquared * 1_eV / (line.GetV0().norm() * p.GetEnergy() * 1_V);
+ geometry::Vector<dimensionless_d> const directionBefore = line.GetV0().normalized();
- geometry::Line lineWithoutB(currentPosition, velocity);
- // determine direction of the particle after adding magnetic field
- // assuming magnetic field does not change during movement
// First Movement
- auto position = currentPosition + velocity * min / 2;
+ // assuming magnetic field does not change during movement
+ auto position = p.GetPosition() + directionBefore * Steplength / 2;
// Change of direction by magnetic field
- geometry::Vector<dimensionless_d> const directionAfter = directionBefore + velocity.cross(magneticfield) *
- min * k;
+ geometry::Vector<dimensionless_d> const directionAfter = directionBefore + directionBefore.cross(magneticfield) *
+ Steplength * k;
// Second Movement
- position = position + directionAfter * velocity.norm() * min / 2;
- if ((position - currentPosition).GetNorm() != 0_m) {
- geometry::Vector<dimensionless_d> const direction = (position - currentPosition).normalized();
- velocity = direction * velocity.norm();
- } // no velocity update for very small steps
- geometry::Line lineWithB(currentPosition, velocity);
-
- return std::make_tuple(geometry::Trajectory<geometry::Line>(lineWithoutB, min),
- geometry::Trajectory<geometry::Line>(lineWithB, min),
- velocity.norm() * min, minIter->second);
+ position = position + directionAfter * Steplength / 2;
+ auto distance = position - p.GetPosition();
+
+ if(distance.norm() == 0_m) {
+ return std::make_tuple(line, position, directionAfter);
+ }
+ geometry::Line newLine(p.GetPosition(), distance.normalized() * line.GetV0().norm());
+ return std::make_tuple(newLine, position, directionAfter);
}
};
--
GitLab