From 7666c11bcf51a23017753aaf34221af3377a5e8d Mon Sep 17 00:00:00 2001
From: Andre Schmidt <schmidt-a@iklx288.ikp.kit.edu>
Date: Wed, 12 Aug 2020 10:15:26 +0200
Subject: [PATCH] update ObservationPlane
---
Framework/Cascade/Cascade.h | 27 ++--
Framework/Geometry/oCoordinateSystem.h | 125 ------------------
.../ObservationPlane/ObservationPlane.cc | 54 ++++++--
Processes/TrackingLine/TrackingLine.h | 50 +++----
4 files changed, 86 insertions(+), 170 deletions(-)
delete mode 100644 Framework/Geometry/oCoordinateSystem.h
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index c7aec5de0..0dceefe91 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -214,16 +214,16 @@ namespace corsika::cascade {
vParticle.GetEnergy() * units::constants::c;
// determine geometric tracking
- auto [step, geomMaxLength, nextVol] = fTracking.GetTrack(vParticle);
+ auto [stepWithoutB, stepWithB, geomMaxLength, nextVol] = fTracking.GetTrack(vParticle);
[[maybe_unused]] auto const& dummy_nextVol = nextVol;
// convert next_step from grammage to length
LengthType const distance_interact =
- currentLogicalNode->GetModelProperties().ArclengthFromGrammage(step,
+ currentLogicalNode->GetModelProperties().ArclengthFromGrammage(stepWithB,
next_interact);
// determine the maximum geometric step length
- LengthType const distance_max = fProcessSequence.MaxStepLength(vParticle, step);
+ LengthType const distance_max = fProcessSequence.MaxStepLength(vParticle, stepWithoutB);
std::cout << "distance_max=" << distance_max << std::endl;
// take minimum of geometry, interaction, decay for next step
@@ -232,23 +232,23 @@ namespace corsika::cascade {
C8LOG_DEBUG("transport particle by : {} m", min_distance / 1_m);
- //determine displacement by the magnetic field
+ // determine displacement by the magnetic field
auto const* currentLogicalVolumeNode = vParticle.GetNode();
int chargeNumber;
- if(corsika::particles::IsNucleus(vParticle.GetPID())) {
+ 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());
+ 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);
// First Movement
- //assuming magnetic field does not change during 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) *
@@ -260,17 +260,20 @@ namespace corsika::cascade {
vParticle.SetMomentum(directionAfter.normalized() * vParticle.GetMomentum().GetNorm());
vParticle.SetPosition(position);
} else {
- vParticle.SetPosition(step.PositionFromArclength(min_distance));
+ vParticle.SetPosition(stepWithoutB.PositionFromArclength(min_distance));
}
// .... also update time, momentum, direction, ...
vParticle.SetTime(vParticle.GetTime() + min_distance / units::constants::c);
std::cout << "New Position: " << vParticle.GetPosition().GetCoordinates() << std::endl;
- step.LimitEndTo(min_distance);
+ // is this necessary?
+ stepWithoutB.LimitEndTo(min_distance);
+ stepWithB.LimitEndTo(min_distance);
// apply all continuous processes on particle + track
- if (process_sequence_.DoContinuous(vParticle, step) ==
- process::EProcessReturn::eParticleAbsorbed) {
+ process::EProcessReturn status = fProcessSequence.DoContinuous(vParticle, stepWithoutB);
+
+ if (status == process::EProcessReturn::eParticleAbsorbed) {
C8LOG_DEBUG("Cascade: delete absorbed particle PID={} E={} GeV",
vParticle.GetPID(), vParticle.GetEnergy() / 1_GeV);
if (!vParticle.isDeleted()) vParticle.Delete();
diff --git a/Framework/Geometry/oCoordinateSystem.h b/Framework/Geometry/oCoordinateSystem.h
deleted file mode 100644
index 614c76da9..000000000
--- a/Framework/Geometry/oCoordinateSystem.h
+++ /dev/null
@@ -1,125 +0,0 @@
-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#ifndef _include_COORDINATESYSTEM_H_
-#define _include_COORDINATESYSTEM_H_
-
-#include <corsika/geometry/QuantityVector.h>
-#include <corsika/units/PhysicalUnits.h>
-#include <corsika/utl/sgn.h>
-#include <Eigen/Dense>
-#include <stdexcept>
-
-typedef Eigen::Transform<double, 3, Eigen::Affine> EigenTransform;
-typedef Eigen::Translation<double, 3> EigenTranslation;
-
-namespace corsika::geometry {
-
- class RootCoordinateSystem;
- template <typename T>
- class Vector;
-
- using corsika::units::si::length_d;
-
- class CoordinateSystem {
- CoordinateSystem const* reference = nullptr;
- EigenTransform transf;
-
- CoordinateSystem(CoordinateSystem const& reference, EigenTransform const& transf)
- : reference(&reference)
- , transf(transf) {}
-
- CoordinateSystem()
- : // for creating the root CS
- transf(EigenTransform::Identity()) {}
-
- protected:
- static auto CreateCS() { return CoordinateSystem(); }
- friend corsika::geometry::RootCoordinateSystem; /// this is the only class that can
- /// create ONE unique root CS
-
- public:
- static EigenTransform GetTransformation(CoordinateSystem const& c1,
- CoordinateSystem const& c2);
-
- auto& operator=(const CoordinateSystem& pCS) {
- reference = pCS.reference;
- transf = pCS.transf;
- return *this;
- }
-
- auto translate(QuantityVector<length_d> vector) const {
- EigenTransform const translation{EigenTranslation(vector.eVector)};
-
- return CoordinateSystem(*this, translation);
- }
-
- template <typename TDim>
- auto RotateToZ(Vector<TDim> vVec) const {
- auto const a = vVec.normalized().GetComponents(*this).eVector;
- auto const a1 = a(0), a2 = a(1);
-
- auto const s = utl::sgn(a(2));
- auto const c = 1 / (1 + s * a(2));
-
- Eigen::Matrix3d A, B;
-
- if (s > 0) {
- A << 1, 0, a1, // comment to prevent clang-format
- 0, 1, a2, // .
- -a1, -a2, 1; // .
- B << -a1 * a1 * c, -a1 * a2 * c, 0, // .
- -a1 * a2 * c, -a2 * a2 * c, 0, // .
- 0, 0, -(a1 * a1 + a2 * a2) * c; // .
-
- } else {
- A << 1, 0, a1, // .
- 0, -1, a2, // .
- a1, -a2, -1; // .
- B << -a1 * a1 * c, +a1 * a2 * c, 0, // .
- -a1 * a2 * c, +a2 * a2 * c, 0, // .
- 0, 0, (a1 * a1 + a2 * a2) * c; // .
- }
-
- return CoordinateSystem(*this, EigenTransform(A + B));
- }
-
- template <typename TDim>
- auto rotate(QuantityVector<TDim> axis, double angle) const {
- if (axis.eVector.isZero()) {
- throw std::runtime_error("null-vector given as axis parameter");
- }
-
- EigenTransform const rotation{Eigen::AngleAxisd(angle, axis.eVector.normalized())};
-
- return CoordinateSystem(*this, rotation);
- }
-
- template <typename TDim>
- auto translateAndRotate(QuantityVector<phys::units::length_d> translation,
- QuantityVector<TDim> axis, double angle) {
- if (axis.eVector.isZero()) {
- throw std::runtime_error("null-vector given as axis parameter");
- }
-
- EigenTransform const transf{Eigen::AngleAxisd(angle, axis.eVector.normalized()) *
- EigenTranslation(translation.eVector)};
-
- return CoordinateSystem(*this, transf);
- }
-
- auto const* GetReference() const { return reference; }
-
- auto const& GetTransform() const { return transf; }
- };
-
-} // namespace corsika::geometry
-
-#endif
diff --git a/Processes/ObservationPlane/ObservationPlane.cc b/Processes/ObservationPlane/ObservationPlane.cc
index 7abff3513..a6f939e33 100644
--- a/Processes/ObservationPlane/ObservationPlane.cc
+++ b/Processes/ObservationPlane/ObservationPlane.cc
@@ -59,20 +59,56 @@ corsika::process::EProcessReturn ObservationPlane::DoContinuous(
}
}
-LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const&,
+LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const& vParticle,
setup::Trajectory const& trajectory) {
- TimeType const timeOfIntersection =
+ int chargeNumber;
+ if (corsika::particles::IsNucleus(vParticle.GetPID())) {
+ chargeNumber = vParticle.GetNuclearZ();
+ } else {
+ chargeNumber = corsika::particles::GetChargeNumber(vParticle.GetPID());
+ }
+ auto const* currentLogicalVolumeNode = vParticle.GetNode();
+ 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) {
+ auto const* currentLogicalVolumeNode = vParticle.GetNode();
+ auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(vParticle.GetPosition());
+ auto k = chargeNumber * corsika::units::constants::cSquared * 1_eV /
+ (velocity.GetSquaredNorm() * vParticle.GetEnergy() * 1_V);
+ LengthType MaxStepLength1 =
+ ( sqrt(velocity.dot(plane_.GetNormal()) * velocity.dot(plane_.GetNormal()) / velocity.GetSquaredNorm() -
+ (plane_.GetNormal().dot(trajectory.GetR0() - plane_.GetCenter()) *
+ plane_.GetNormal().dot(velocity.cross(magneticfield)) * 2 * k)) -
+ velocity.dot(plane_.GetNormal()) / velocity.GetNorm() ) /
+ (plane_.GetNormal().dot(velocity.cross(magneticfield)) * k);
+ LengthType MaxStepLength2 =
+ ( - sqrt(velocity.dot(plane_.GetNormal()) * velocity.dot(plane_.GetNormal()) / velocity.GetSquaredNorm() -
+ (plane_.GetNormal().dot(trajectory.GetR0() - plane_.GetCenter()) *
+ plane_.GetNormal().dot(velocity.cross(magneticfield)) * 2 * k)) -
+ velocity.dot(plane_.GetNormal()) / velocity.GetNorm() ) /
+ (plane_.GetNormal().dot(velocity.cross(magneticfield)) * k);
+ std::cout << "Test: " << MaxStepLength1 << " " << MaxStepLength2 << std::endl;
+ if (MaxStepLength1 < 0_m && MaxStepLength2 < 0_m) {
+ return std::numeric_limits<double>::infinity() * 1_m;
+ } else if (MaxStepLength1 < 0_m || MaxStepLength2 < MaxStepLength1) {
+ return MaxStepLength2;
+ } else if (MaxStepLength2 < 0_m || MaxStepLength1 < MaxStepLength2) {
+ return MaxStepLength1;
+ }
+ } else {
+ TimeType const timeOfIntersection =
(plane_.GetCenter() - trajectory.GetR0()).dot(plane_.GetNormal()) /
trajectory.GetV0().dot(plane_.GetNormal());
- if (timeOfIntersection < TimeType::zero()) {
- return std::numeric_limits<double>::infinity() * 1_m;
- }
+ if (timeOfIntersection < TimeType::zero()) {
+ return std::numeric_limits<double>::infinity() * 1_m;
+ }
- auto const pointOfIntersection = trajectory.GetPosition(timeOfIntersection);
- auto dist = (trajectory.GetR0() - pointOfIntersection).norm() * 1.0001;
- C8LOG_TRACE("ObservationPlane::MaxStepLength l={} m", dist / 1_m);
- return dist;
+ auto const pointOfIntersection = trajectory.GetPosition(timeOfIntersection);
+ return (trajectory.GetR0() - pointOfIntersection).norm() * 1.0001;
+ //why is it *1.0001? should i do that too?
+ }
}
void ObservationPlane::ShowResults() const {
diff --git a/Processes/TrackingLine/TrackingLine.h b/Processes/TrackingLine/TrackingLine.h
index 138c2f191..c16dd44b1 100644
--- a/Processes/TrackingLine/TrackingLine.h
+++ b/Processes/TrackingLine/TrackingLine.h
@@ -77,7 +77,7 @@ namespace corsika::process {
//charge of the particle
int chargeNumber;
- if(corsika::particles::IsNucleus(p.GetPID())) {
+ if (corsika::particles::IsNucleus(p.GetPID())) {
chargeNumber = p.GetNuclearZ();
} else {
chargeNumber = corsika::particles::GetChargeNumber(p.GetPID());
@@ -96,9 +96,9 @@ namespace corsika::process {
volume); // for the moment we are a bit bold here and assume
// everything is a sphere, crashes with exception if not
- //creating Line with magnetic field
- if(chargeNumber != 0) {
- //determine steplength to next volume
+ // creating Line with magnetic field
+ if (chargeNumber != 0) {
+ // determine steplength to next volume
double a = ((directionBefore.cross(magneticfield)).dot(currentPosition - sphere.GetCenter()) * k + 1) * 4 /
(1_m * 1_m * (directionBefore.cross(magneticfield)).GetSquaredNorm() * k * k);
double b = directionBefore.dot(currentPosition - sphere.GetCenter()) * 8 /
@@ -108,22 +108,22 @@ namespace corsika::process {
((directionBefore.cross(magneticfield)).GetSquaredNorm() * k * k * 1_m * 1_m * 1_m * 1_m);
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) {
+ for (int i = 0; i < 4; i++) {
+ if (solutions[i].imag() == 0 && solutions[i].real() > 0) {
tmp.push_back(solutions[i].real());
}
}
LengthType Steplength;
- if(tmp.size() > 0) {
+ if (tmp.size() > 0) {
Steplength = 1_m * *std::min_element(tmp.begin(),tmp.end());
std::cout << "s = " << Steplength << std::endl;
} else {
std::cout << "no intersection (1)!" << std::endl;
- //what to do when this happens?
+ // what to do when this happens? (very unlikely)
}
// First Movement
- //assuming magnetic field does not change during 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) *
@@ -133,8 +133,8 @@ namespace corsika::process {
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
- //line has some errors
+ } // instead of changing the line with magnetic field, the TimeOfIntersection() could be changed
+ // line has some errors
geometry::Line line(currentPosition, velocity1);
if (auto opt = TimeOfIntersection(line, sphere); opt.has_value()) {
@@ -156,9 +156,9 @@ namespace corsika::process {
// for the moment we are a bit bold here and assume
// everything is a sphere, crashes with exception if not
- //creating Line with magnetic field
- if(chargeNumber != 0) {
- //determine steplength to next volume
+ // creating Line with magnetic field
+ if (chargeNumber != 0) {
+ // determine steplength to next volume
double a = ((directionBefore.cross(magneticfield)).dot(currentPosition - sphere.GetCenter()) * k + 1) * 4 /
(1_m * 1_m * (directionBefore.cross(magneticfield)).GetSquaredNorm() * k * k);
double b = directionBefore.dot(currentPosition - sphere.GetCenter()) * 8 /
@@ -168,22 +168,22 @@ namespace corsika::process {
((directionBefore.cross(magneticfield)).GetSquaredNorm() * k * k * 1_m * 1_m * 1_m * 1_m);
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) {
+ for (int i = 0; i < 4; i++) {
+ if (solutions[i].imag() == 0 && solutions[i].real() > 0) {
tmp.push_back(solutions[i].real());
}
}
LengthType Steplength;
- if(tmp.size() > 0) {
+ if (tmp.size() > 0) {
Steplength = 1_m * *std::min_element(tmp.begin(),tmp.end());
std::cout << "s = " << Steplength << std::endl;
} else {
std::cout << "no intersection (2)!" << std::endl;
- //what to do when this happens?
+ // what to do when this happens? (very unlikely)
}
// First Movement
- //assuming magnetic field does not change during 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) *
@@ -193,7 +193,7 @@ namespace corsika::process {
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
+ } // instead of changing the line with magnetic field, the TimeOfIntersection() could be changed
geometry::Line line(currentPosition, velocity2);
[[maybe_unused]] auto const [t1, t2] = *TimeOfIntersection(line, sphere);
@@ -219,9 +219,10 @@ namespace corsika::process {
<< min
// << " " << minIter->second->GetModelProperties().GetName()
<< std::endl;
-
+
+ geometry::Line lineWithoutB(currentPosition, velocity);
// determine direction of the particle after adding magnetic field
- //assuming magnetic field does not change during movement
+ // assuming magnetic field does not change during movement
// First Movement
auto position = currentPosition + velocity * min / 2;
// Change of direction by magnetic field
@@ -232,9 +233,10 @@ namespace corsika::process {
geometry::Vector<dimensionless_d> const direction = (position - currentPosition) /
(position - currentPosition).GetNorm();
velocity = direction * velocity.GetNorm();
- geometry::Line line(currentPosition, velocity);
+ geometry::Line lineWithB(currentPosition, velocity);
- return std::make_tuple(geometry::Trajectory<geometry::Line>(line, min),
+ return std::make_tuple(geometry::Trajectory<geometry::Line>(lineWithoutB, min),
+ geometry::Trajectory<geometry::Line>(lineWithB, min),
velocity.norm() * min, minIter->second);
}
};
--
GitLab