IAP GITLAB

Commit cf2d0c5f authored by Ralf Ulrich's avatar Ralf Ulrich

fixed missing curved-plane intersection code and a few typos

parent c36652aa
...@@ -18,9 +18,6 @@ ...@@ -18,9 +18,6 @@
#include <corsika/framework/stack/SecondaryView.hpp> #include <corsika/framework/stack/SecondaryView.hpp>
#include <corsika/media/Environment.hpp> #include <corsika/media/Environment.hpp>
#include <corsika/setup/SetupStack.hpp>
#include <corsika/setup/SetupTrajectory.hpp>
#include <cassert> #include <cassert>
#include <cmath> #include <cmath>
#include <iostream> #include <iostream>
...@@ -139,7 +136,6 @@ namespace corsika { ...@@ -139,7 +136,6 @@ namespace corsika {
distance_interact / 1_m, continuous_max_dist / 1_m); distance_interact / 1_m, continuous_max_dist / 1_m);
// here the particle is actually moved along the trajectory to new position: // here the particle is actually moved along the trajectory to new position:
// std::visit(setup::ParticleUpdate<particle_type>{vParticle}, step);
step.setLength(min_distance); step.setLength(min_distance);
vParticle.setPosition(step.getPosition(1)); vParticle.setPosition(step.getPosition(1));
// assumption: tracking does not change absolute momentum: // assumption: tracking does not change absolute momentum:
......
...@@ -30,7 +30,7 @@ namespace corsika { ...@@ -30,7 +30,7 @@ namespace corsika {
template <typename TParticle> template <typename TParticle>
inline auto make_LeapFrogStep(TParticle const& particle, LengthType steplength) { inline auto make_LeapFrogStep(TParticle const& particle, LengthType steplength) {
if (particle.getMomentum().norm() == 0_GeV) { if (particle.getMomentum().getNorm() == 0_GeV) {
return std::make_tuple(particle.getPosition(), particle.getMomentum() / 1_GeV, return std::make_tuple(particle.getPosition(), particle.getMomentum() / 1_GeV,
double(0)); double(0));
} // charge of the particle } // charge of the particle
...@@ -225,6 +225,81 @@ namespace corsika { ...@@ -225,6 +225,81 @@ namespace corsika {
return Intersections(d_enter / absVelocity, d_exit / absVelocity); return Intersections(d_enter / absVelocity, d_exit / absVelocity);
} }
template <typename TParticle, typename TMedium>
Intersections Tracking::intersect(TParticle const& particle, Plane const& plane,
TMedium const& medium) {
int chargeNumber;
if (is_nucleus(particle.getPID())) {
chargeNumber = particle.getNuclearZ();
} else {
chargeNumber = get_charge_number(particle.getPID());
}
auto const* currentLogicalVolumeNode = particle.getNode();
auto magneticfield = currentLogicalVolumeNode->getModelProperties().getMagneticField(
particle.getPosition());
VelocityVector const velocity =
particle.getMomentum() / particle.getEnergy() * constants::c;
auto const absVelocity = velocity.getNorm();
DirectionVector const direction =
velocity.normalized(); // determine steplength to next volume
Point const position = particle.getPosition();
if (chargeNumber != 0 &&
abs(plane.getNormal().dot(velocity.cross(magneticfield))) *
1_s / 1_m / 1_T > 1e-6) {
auto const* currentLogicalVolumeNode = particle.getNode();
auto magneticfield = currentLogicalVolumeNode->getModelProperties().getMagneticField(
particle.getPosition());
auto k = chargeNumber * constants::c * 1_eV /
(particle.getMomentum().getNorm() * 1_V);
if (direction.dot(plane.getNormal()) * direction.dot(plane.getNormal()) -
(plane.getNormal().dot(position - plane.getCenter()) *
plane.getNormal().dot(direction.cross(magneticfield)) * 2 * k) <
0) {
return Intersections(std::numeric_limits<double>::infinity() * 1_s);
}
LengthType const MaxStepLength1 =
(sqrt(direction.dot(plane.getNormal()) * direction.dot(plane.getNormal()) -
(plane.getNormal().dot(position - plane.getCenter()) *
plane.getNormal().dot(direction.cross(magneticfield)) * 2 * k)) -
direction.dot(plane.getNormal()) / direction.getNorm()) /
(plane.getNormal().dot(direction.cross(magneticfield)) * k);
LengthType const MaxStepLength2 =
(-sqrt(
direction.dot(plane.getNormal()) * direction.dot(plane.getNormal()) -
(plane.getNormal().dot(position - plane.getCenter()) *
plane.getNormal().dot(direction.cross(magneticfield)) * 2 * k)) -
direction.dot(plane.getNormal()) / direction.getNorm()) /
(plane.getNormal().dot(direction.cross(magneticfield)) * k);
if (MaxStepLength1 <= 0_m && MaxStepLength2 <= 0_m) {
return Intersections(std::numeric_limits<double>::infinity() * 1_s);
} else if (MaxStepLength1 <= 0_m || MaxStepLength2 < MaxStepLength1) {
CORSIKA_LOG_TRACE(" steplength to obs plane 2: {} ", MaxStepLength2);
return Intersections(MaxStepLength2 *
(direction + direction.cross(magneticfield) * MaxStepLength2 * k / 2)
.getNorm() / absVelocity);
} else if (MaxStepLength2 <= 0_m || MaxStepLength1 < MaxStepLength2) {
CORSIKA_LOG_TRACE(" steplength to obs plane 2: {} ", MaxStepLength1);
return Intersections(MaxStepLength1 *
(direction + direction.cross(magneticfield) * MaxStepLength2 * k / 2)
.getNorm() / absVelocity);
}
CORSIKA_LOG_WARN("Particle wasn't tracked with curved trajectory -> straight");
} // end if curved-tracking
return tracking_line::Tracking::intersect(particle, plane, medium);
}
template <typename TParticle, typename TBaseNodeType> template <typename TParticle, typename TBaseNodeType>
inline Intersections Tracking::intersect(const TParticle& particle, inline Intersections Tracking::intersect(const TParticle& particle,
const TBaseNodeType& volumeNode) { const TBaseNodeType& volumeNode) {
......
...@@ -9,8 +9,6 @@ ...@@ -9,8 +9,6 @@
#pragma once #pragma once
#include <corsika/framework/geometry/Line.hpp> #include <corsika/framework/geometry/Line.hpp>
#include <corsika/framework/geometry/Plane.hpp>
#include <corsika/framework/geometry/Sphere.hpp>
#include <corsika/framework/geometry/StraightTrajectory.hpp> #include <corsika/framework/geometry/StraightTrajectory.hpp>
#include <corsika/framework/geometry/Vector.hpp> #include <corsika/framework/geometry/Vector.hpp>
#include <corsika/framework/core/ParticleProperties.hpp> #include <corsika/framework/core/ParticleProperties.hpp>
......
...@@ -61,12 +61,16 @@ namespace corsika { ...@@ -61,12 +61,16 @@ namespace corsika {
auto getTrack(TParticle const& particle); auto getTrack(TParticle const& particle);
template <typename TParticle, typename TMedium> template <typename TParticle, typename TMedium>
static Intersections intersect(const TParticle& particle, const Sphere& sphere, static Intersections intersect(TParticle const& particle, Sphere const& sphere,
const TMedium& medium); TMedium const& medium);
template <typename TParticle, typename TBaseNodeType> template <typename TParticle, typename TBaseNodeType>
static Intersections intersect(const TParticle& particle, static Intersections intersect(TParticle const& particle,
const TBaseNodeType& volumeNode); TBaseNodeType const& volumeNode);
template <typename TParticle, typename TMedium>
static Intersections intersect(TParticle const& particle, Plane const& plane,
TMedium const& medium);
protected: protected:
/** /**
......
...@@ -38,8 +38,8 @@ namespace corsika::setup { ...@@ -38,8 +38,8 @@ namespace corsika::setup {
The default tracking algorithm. The default tracking algorithm.
*/ */
typedef corsika::process::tracking_leapfrog_curved::Tracking Tracking; typedef corsika::tracking_leapfrog_curved::Tracking Tracking;
// typedef corsika::process::tracking_leapfrog_straight::Tracking Tracking; // typedef corsika::tracking_leapfrog_straight::Tracking Tracking;
// typedef corsika::tracking_line::Tracking Tracking; // typedef corsika::tracking_line::Tracking Tracking;
/** /**
...@@ -47,6 +47,6 @@ namespace corsika::setup { ...@@ -47,6 +47,6 @@ namespace corsika::setup {
*/ */
/// definition of Trajectory base class, to be used in tracking and cascades /// definition of Trajectory base class, to be used in tracking and cascades
// typedef StraightTrajectory Trajectory; // typedef StraightTrajectory Trajectory;
typedef corsika::geometry::LeapFrogTrajectory Trajectory; typedef corsika::LeapFrogTrajectory Trajectory;
} // namespace corsika::setup } // namespace corsika::setup
...@@ -13,33 +13,33 @@ ...@@ -13,33 +13,33 @@
#include <corsika/framework/geometry/Helix.hpp> #include <corsika/framework/geometry/Helix.hpp>
#include <corsika/framework/geometry/StraightTrajectory.hpp> #include <corsika/framework/geometry/StraightTrajectory.hpp>
// #include <corsika/framework/geometry/LeapFrogTrajectory.hpp> #include <corsika/framework/geometry/LeapFrogTrajectory.hpp>
// #include <corsika/modules/TrackingLine.hpp> //#include <corsika/modules/TrackingLine.hpp>
// #include <corsika/modules/TrackingCurved.hpp> // simple leap-frog implementation //#include <corsika/modules/TrackingCurved.hpp> // simple leap-frog implementation
// #include <corsika/modules/TrackingLeapFrog.hpp> // more complete leap-frog //#include <corsika/modules/TrackingLeapFrog.hpp> // more complete leap-frog
// implementation // implementation
namespace corsika::setup::testing { namespace corsika::setup::testing {
template <typename TTrack> template <typename TTrack>
TTrack make_track(Line const& line, TimeType const tEnd); TTrack make_track(Line const line,
TimeType const tEnd = std::numeric_limits<TimeType::value_type>::infinity() * 1_s);
template <> template <>
inline StraightTrajectory make_track<StraightTrajectory>(Line const& line, inline StraightTrajectory make_track<StraightTrajectory>(Line const line,
TimeType const tEnd) { TimeType const tEnd) {
return StraightTrajectory(line, tEnd); return StraightTrajectory(line, tEnd);
} }
/* template <>
template <> inline LeapFrogTrajectory make_track<LeapFrogTrajectory>(Line const line,
inline LeapFrogTrajectory make_track<LeapFrogTrajectory>(Line const& line, TimeType const tEnd) {
TimeType const tEnd) {
auto const k = square(0_m) / (square(1_s) * 1_V);
auto const k = square(0_m) / (square(1_s) * 1_V); return LeapFrogTrajectory(
return LeapFrogTrajectory( line.getStartPoint(), line.getVelocity(),
line.getStartPoint(), line.getVelocity(), MagneticFieldVector{line.getStartPoint().getCoordinateSystem(), 0_T, 0_T, 0_T},
MagneticFieldVector{line.getStartPoint().getCoordinateSystem(), 0_T, 0_T, 0_T}, k, tEnd);
k, tEnd);
} }
*/
} // namespace corsika::setup::testing } // namespace corsika::setup::testing
...@@ -24,6 +24,8 @@ ...@@ -24,6 +24,8 @@
#include <corsika/media/HomogeneousMedium.hpp> #include <corsika/media/HomogeneousMedium.hpp>
#include <corsika/media/NuclearComposition.hpp> #include <corsika/media/NuclearComposition.hpp>
#include <SetupTestTrajectory.hpp>
#include <catch2/catch.hpp> #include <catch2/catch.hpp>
using namespace corsika; using namespace corsika;
...@@ -68,14 +70,11 @@ public: ...@@ -68,14 +70,11 @@ public:
auto getTrack(TParticle const& particle) { auto getTrack(TParticle const& particle) {
VelocityVector const initialVelocity = VelocityVector const initialVelocity =
particle.getMomentum() / particle.getEnergy() * constants::c; particle.getMomentum() / particle.getEnergy() * constants::c;
Line const theLine = Line(particle.getPosition(), initialVelocity);
TimeType const tEnd = std::numeric_limits<TimeType::value_type>::infinity() * 1_s;
return std::make_tuple( return std::make_tuple(
StraightTrajectory( corsika::setup::testing::make_track<setup::Trajectory>(theLine, tEnd),
Line(particle.getPosition(), initialVelocity), // trajectory: just go ahead forever
std::numeric_limits<TimeType::value_type>::infinity() * 1_s), // trajectory,
// just
// go
// ahead
// forever
particle.getNode()); // next volume node particle.getNode()); // next volume node
} }
}; };
......
...@@ -101,7 +101,7 @@ TEST_CASE("MediumPropertyModel w/ Homogeneous") { ...@@ -101,7 +101,7 @@ TEST_CASE("MediumPropertyModel w/ Homogeneous") {
auto const tEnd = 1_s; auto const tEnd = 1_s;
// and the associated trajectory // and the associated trajectory
setup::Trajectory const trajectory(line, tEnd); setup::Trajectory const trajectory = corsika::setup::testing::make_track<setup::Trajectory>(line, tEnd);
// and check the integrated grammage // and check the integrated grammage
CHECK((medium.getIntegratedGrammage(trajectory, 3_m) / (density * 3_m)) == Approx(1)); CHECK((medium.getIntegratedGrammage(trajectory, 3_m) / (density * 3_m)) == Approx(1));
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment