From 9bc74cc694b46a628ba21811f9e5b22ce5ff7200 Mon Sep 17 00:00:00 2001
From: Andre Schmidt <schmidt-a@iklx288.ikp.kit.edu>
Date: Thu, 6 Aug 2020 11:02:14 +0200
Subject: [PATCH] change Steplength
---
Documentation/Examples/vertical_EAS.cc | 2 +
Framework/Cascade/Cascade.h | 40 ++++++++++--
Framework/Utilities/CMakeLists.txt | 2 +
Processes/TrackingLine/TrackingLine.h | 74 ++++++++++++----------
Processes/TrackingLine/testTrackingLine.cc | 2 +-
5 files changed, 81 insertions(+), 39 deletions(-)
diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc
index 0f4d57c87..25f246748 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -34,6 +34,7 @@
#include <corsika/process/particle_cut/ParticleCut.h>
#include <corsika/process/proposal/ContinuousProcess.h>
#include <corsika/process/proposal/Interaction.h>
+#include <corsika/process/track_writer/TrackWriter.h>
#include <corsika/process/pythia/Decay.h>
#include <corsika/process/sibyll/Decay.h>
#include <corsika/process/sibyll/Interaction.h>
@@ -226,6 +227,7 @@ int main(int argc, char** argv) {
process::on_shell_check::OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);
+ process::track_writer::TrackWriter trackWriter("tracks.dat");
process::longitudinal_profile::LongitudinalProfile longprof{showerAxis};
Plane const obsPlane(showerCore, Vector<dimensionless_d>(rootCS, {0., 0., 1.}));
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index 4481e8a24..fe5ecaaff 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 [step, geomMaxLength, nextVol, magMaxLength, directionBefore, directionAfter] = fTracking.GetTrack(vParticle);
+ auto [step, geomMaxLength, nextVol] = fTracking.GetTrack(vParticle);
[[maybe_unused]] auto const& dummy_nextVol = nextVol;
// convert next_step from grammage to length
@@ -226,18 +226,46 @@ namespace corsika::cascade {
LengthType const distance_max = fProcessSequence.MaxStepLength(vParticle, step);
std::cout << "distance_max=" << distance_max << std::endl;
- // take minimum of geometry, interaction, decay for next step
+ // take minimum of geometry, interaction, decay, magnetic field for next step
auto const min_distance = std::min(
- {distance_interact, distance_decay, distance_max, geomMaxLength, magMaxLength});
+ {distance_interact, distance_decay, distance_max, geomMaxLength});
C8LOG_DEBUG("transport particle by : {} m", min_distance / 1_m);
+ //determine displacement by the magnetic field
+ auto const* currentLogicalVolumeNode = vParticle.GetNode();
+ int chargeNumber;
+ if(corsika::particles::IsNucleus(vParticle.GetPID())) {
+ chargeNumber = vParticle.GetNuclearZ();
+ } else {
+ chargeNumber = corsika::particles::GetChargeNumber(vParticle.GetPID());
+ }
+ if(chargeNumber != 0) {
+ 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();
+ auto k = chargeNumber * corsika::units::constants::cSquared * 1_eV /
+ (velocity.GetNorm() * vParticle.GetEnergy() * 1_V);
+ LengthType const steplength = min_distance /
+ (directionBefore + directionBefore.cross(magneticfield) * k / 2).GetNorm();
+ // First Movement
+ //assuming magnetic field does not change during movement
+ auto position = vParticle.GetPosition() + 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 - vParticle.GetPosition()) /
+ (position - vParticle.GetPosition()).GetNorm();
+ vParticle.SetMomentum(direction * vParticle.GetMomentum().GetNorm());
+ }
+
// here the particle is actually moved along the trajectory to new position:
// std::visit(setup::ParticleUpdate<Particle>{vParticle}, step);
vParticle.SetPosition(step.PositionFromArclength(min_distance));
// .... also update time, momentum, direction, ...
- vParticle.SetMomentum((directionBefore * (1 - min_distance / magMaxLength) +
- directionAfter * min_distance /magMaxLength) * vParticle.GetMomentum().GetNorm());
vParticle.SetTime(vParticle.GetTime() + min_distance / units::constants::c);
step.LimitEndTo(min_distance);
@@ -263,7 +291,7 @@ namespace corsika::cascade {
TStackView secondaries(vParticle);
- if (min_distance != distance_max && min_distance != magMaxLength) {
+ if (min_distance != distance_max) {
/*
Create SecondaryView object on Stack. The data container
remains untouched and identical, and 'projectil' is identical
diff --git a/Framework/Utilities/CMakeLists.txt b/Framework/Utilities/CMakeLists.txt
index 95170d4a3..3373c0171 100644
--- a/Framework/Utilities/CMakeLists.txt
+++ b/Framework/Utilities/CMakeLists.txt
@@ -20,6 +20,7 @@ set (
UTILITIES_SOURCES
COMBoost.cc
CorsikaData.cc
+ quartic.cpp
${CORSIKA_FENV})
set (
@@ -31,6 +32,7 @@ set (
sgn.h
CorsikaFenv.h
MetaProgramming.h
+ quartic.h
)
set (
diff --git a/Processes/TrackingLine/TrackingLine.h b/Processes/TrackingLine/TrackingLine.h
index 5282f5a57..de024f564 100644
--- a/Processes/TrackingLine/TrackingLine.h
+++ b/Processes/TrackingLine/TrackingLine.h
@@ -15,8 +15,7 @@
#include <corsika/geometry/Vector.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/units/PhysicalUnits.h>
-#include <corsika/logging/Logging.h>
-
+#include <corsika/utl/quartic.h>
#include <optional>
#include <type_traits>
#include <utility>
@@ -60,7 +59,7 @@ namespace corsika::process {
<< std::endl;
std::cout << "TrackingLine E: " << p.GetEnergy() / 1_GeV << " GeV" << std::endl;
- // determine velocity after adding magnetic field
+ // determine direction of the particle after adding magnetic field
auto const* currentLogicalVolumeNode = p.GetNode();
int chargeNumber;
if(corsika::particles::IsNucleus(p.GetPID())) {
@@ -68,39 +67,51 @@ namespace corsika::process {
} else {
chargeNumber = corsika::particles::GetChargeNumber(p.GetPID());
}
- geometry::Vector<dimensionless_d> const directionBefore = velocity.normalized();
- auto magMaxLength = 1_m/0;
- auto directionAfter = directionBefore;
if(chargeNumber != 0) {
auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(currentPosition);
std::cout << "TrackingLine B: " << magneticfield.GetComponents() / 1_uT << " uT " << std::endl;
- geometry::Vector<SpeedType::dimension_type> const velocityVerticalMag = velocity -
- velocity.parallelProjectionOnto(magneticfield);
- LengthType const gyroradius = p.GetEnergy() * velocityVerticalMag.GetNorm() * 1_V /
- (corsika::units::constants::cSquared * abs(chargeNumber) *
- magneticfield.GetNorm() * 1_eV);
- //steplength depending on how exact it should be
- LengthType const Steplength = 0.01 * gyroradius;
- // First Movement
- auto position = currentPosition + directionBefore * Steplength / 2;
- // Change of direction by magnetic field at position
- magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(position);
- directionAfter = directionBefore + directionBefore.cross(magneticfield) * chargeNumber *
- Steplength * corsika::units::constants::cSquared * 1_eV /
- (p.GetEnergy() * velocity.GetNorm() * 1_V);
- // Second Movement
- position = position + directionAfter * Steplength / 2;
- magMaxLength = (position - currentPosition).GetNorm();
- geometry::Vector<dimensionless_d> const direction = (position - currentPosition) /
- magMaxLength;
- velocity = direction * velocity.GetNorm();
- std::cout << "TrackingLine p: " << (direction * p.GetMomentum().GetNorm()).GetComponents() / 1_GeV
- << " GeV " << 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();
+ //determine steplength to next volume
+ double a = ((directionBefore.cross(magneticfield)).dot(currentPosition - .GetCenter()) * k + 1) /
+ ((directionBefore.cross(magneticfield)).GetSquaredNorm() * k^2 * 1_m * 1_m / 4);
+ double b = directionBefore.dot(currentPosition - .GetCenter()) * 2 /
+ ((directionBefore.cross(magneticfield)).GetSquaredNorm() * k^2 * 1_m / 4);
+ double c = ((currentPosition - .GetCenter()).GetSquaredNorm() - .GetRadius^2) /
+ ((directionBefore.cross(magneticfield)).GetSquaredNorm() * k^2 / 4);
+ std::complex<double>* solutions = solve_quartic(0, a, b, c);
+ std::vector<double> tmp;
+ for(int i = 0; i < 4; i++) {
+ if(solutions[i].imag() == 0 && solutions[i].real() > 0) {
+ tmp.push_back(solutions[i].real());
+ }
+ }
+ LengthType const Steplength;
+ if(tmp.size() > 0) {
+ Steplength = *std::min_element(tmp.begin(),tmp.end()) * 1_m;
+ std::cout << "s = " << Steplength << std::endl;
+ } else {
+ std::cout << "no intersection with anything!" << std::endl;
+ //what to do when this happens?
+ }
+
+ // 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();
+ velocity = direction * velocity.GetNorm();
+ std::cout << "TrackingLine p: " << (direction * p.GetMomentum().GetNorm()).GetComponents() / 1_GeV
+ << " GeV " << std::endl;
} else {
std::cout << "TrackingLine p: " << p.GetMomentum().GetComponents() / 1_GeV
- << " GeV " << std::endl;
+ << " GeV " << std::endl;
}
-
std::cout << "TrackingLine v: " << velocity.GetComponents() << std::endl;
geometry::Line line(currentPosition, velocity);
@@ -166,8 +177,7 @@ namespace corsika::process {
C8LOG_DEBUG(" t-intersect: {} ", min);
return std::make_tuple(geometry::Trajectory<geometry::Line>(line, min),
- velocity.norm() * min, minIter->second, magMaxLength,
- directionBefore, directionAfter);
+ velocity.norm() * min, minIter->second);
}
};
diff --git a/Processes/TrackingLine/testTrackingLine.cc b/Processes/TrackingLine/testTrackingLine.cc
index 4cfa2bd60..a72e870b7 100644
--- a/Processes/TrackingLine/testTrackingLine.cc
+++ b/Processes/TrackingLine/testTrackingLine.cc
@@ -86,7 +86,7 @@ TEST_CASE("TrackingLine") {
0_m / second, 1_m / second);
Line line(origin, v);
- auto const [traj, geomMaxLength, nextVol] = tracking.GetTrack(p);
+ auto const [traj, geomMaxLength, nextVol, magMaxLength, beforeDirection, afterDirection] = tracking.GetTrack(p);
[[maybe_unused]] auto dummy_geomMaxLength = geomMaxLength;
[[maybe_unused]] auto dummy_nextVol = nextVol;
--
GitLab