From 17cee635dbc9ffe3f00ecb13f77da08e1510afbb Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@kit.edu>
Date: Thu, 5 Mar 2020 17:38:26 +0100
Subject: [PATCH] "transparent" observation plane
---
Processes/ObservationPlane/CMakeLists.txt | 7 ++-
.../ObservationPlane/ObservationPlane.cc | 46 ++++++++++++-------
Processes/ObservationPlane/ObservationPlane.h | 7 +--
.../ObservationPlane/testObservationPlane.cc | 28 +++++------
4 files changed, 51 insertions(+), 37 deletions(-)
diff --git a/Processes/ObservationPlane/CMakeLists.txt b/Processes/ObservationPlane/CMakeLists.txt
index 2707887a1..e1c3c4179 100644
--- a/Processes/ObservationPlane/CMakeLists.txt
+++ b/Processes/ObservationPlane/CMakeLists.txt
@@ -16,7 +16,7 @@ set (
add_library (ProcessObservationPlane STATIC ${MODEL_SOURCES})
CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessObservationPlane ${MODEL_NAMESPACE} ${MODEL_HEADERS})
-# target dependencies on other libraries (also the header onlys)
+#target dependencies on other libraries(also the header onlys)
target_link_libraries (
ProcessObservationPlane
CORSIKAgeometry
@@ -32,8 +32,8 @@ target_include_directories (
install (FILES ${MODEL_HEADERS} DESTINATION include/${MODEL_NAMESPACE})
-# --------------------
-# code unit testing
+#-- -- -- -- -- -- -- -- -- --
+#code unit testing
CORSIKA_ADD_TEST(testObservationPlane)
target_link_libraries (
testObservationPlane
@@ -41,4 +41,3 @@ target_link_libraries (
CORSIKAstackinterface
CORSIKAthirdparty # for catch2
)
-
diff --git a/Processes/ObservationPlane/ObservationPlane.cc b/Processes/ObservationPlane/ObservationPlane.cc
index 69c2f1382..8269294cc 100644
--- a/Processes/ObservationPlane/ObservationPlane.cc
+++ b/Processes/ObservationPlane/ObservationPlane.cc
@@ -15,36 +15,48 @@
using namespace corsika::process::observation_plane;
using namespace corsika::units::si;
-ObservationPlane::ObservationPlane(geometry::Plane const& vObsPlane,
- std::string const& vFilename)
- : fObsPlane(vObsPlane)
- , fOutputStream(vFilename) {
- fOutputStream << "#PDG code, energy / eV, distance to center / m" << std::endl;
+ObservationPlane::ObservationPlane(geometry::Plane const& obsPlane,
+ std::string const& filename, bool deleteOnHit)
+ : plane_(obsPlane)
+ , outputStream_(filename)
+ , deleteOnHit_(deleteOnHit) {
+ outputStream_ << "#PDG code, energy / eV, distance to center / m" << std::endl;
}
corsika::process::EProcessReturn ObservationPlane::DoContinuous(
- setup::Stack::ParticleType const& vParticle, setup::Trajectory const& vTrajectory) {
+ setup::Stack::ParticleType const& particle, setup::Trajectory const& trajectory) {
+ TimeType const timeOfIntersection =
+ (plane_.GetCenter() - trajectory.GetR0()).dot(plane_.GetNormal()) /
+ trajectory.GetV0().dot(plane_.GetNormal());
- if (fObsPlane.IsAbove(vTrajectory.GetPosition(1.0001))) {
+ if (timeOfIntersection < TimeType::zero()) { return process::EProcessReturn::eOk; }
+
+ if (plane_.IsAbove(trajectory.GetR0()) == plane_.IsAbove(trajectory.GetPosition(1))) {
return process::EProcessReturn::eOk;
}
- fOutputStream << static_cast<int>(particles::GetPDG(vParticle.GetPID())) << ' '
- << vParticle.GetEnergy() * (1 / 1_eV) << ' '
- << (vTrajectory.GetPosition(1) - fObsPlane.GetCenter()).norm() / 1_m
+ outputStream_ << static_cast<int>(particles::GetPDG(particle.GetPID())) << ' '
+ << particle.GetEnergy() * (1 / 1_eV) << ' '
+ << (trajectory.GetPosition(1) - plane_.GetCenter()).norm() / 1_m
<< std::endl;
- return process::EProcessReturn::eParticleAbsorbed;
+ if (deleteOnHit_) {
+ return process::EProcessReturn::eParticleAbsorbed;
+ } else {
+ return process::EProcessReturn::eOk;
+ }
}
LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const&,
- setup::Trajectory const& vTrajectory) {
- if (!fObsPlane.IsAbove(vTrajectory.GetR0())) {
+ setup::Trajectory const& trajectory) {
+ 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;
}
- auto const pointOfIntersection = vTrajectory.GetPosition(
- (fObsPlane.GetCenter() - vTrajectory.GetR0()).dot(fObsPlane.GetNormal()) /
- vTrajectory.GetV0().dot(fObsPlane.GetNormal()));
- return (vTrajectory.GetR0() - pointOfIntersection).norm();
+ auto const pointOfIntersection = trajectory.GetPosition(timeOfIntersection);
+ return (trajectory.GetR0() - pointOfIntersection).norm() * 1.0001;
}
diff --git a/Processes/ObservationPlane/ObservationPlane.h b/Processes/ObservationPlane/ObservationPlane.h
index 9ecec9f63..900b7a589 100644
--- a/Processes/ObservationPlane/ObservationPlane.h
+++ b/Processes/ObservationPlane/ObservationPlane.h
@@ -29,7 +29,7 @@ namespace corsika::process::observation_plane {
class ObservationPlane : public corsika::process::ContinuousProcess<ObservationPlane> {
public:
- ObservationPlane(geometry::Plane const& vObsPlane, std::string const& vFilename);
+ ObservationPlane(geometry::Plane const&, std::string const&, bool = true);
void Init() {}
corsika::process::EProcessReturn DoContinuous(
@@ -41,8 +41,9 @@ namespace corsika::process::observation_plane {
corsika::setup::Trajectory const& vTrajectory);
private:
- geometry::Plane const fObsPlane;
- std::ofstream fOutputStream;
+ geometry::Plane const plane_;
+ std::ofstream outputStream_;
+ bool const deleteOnHit_;
};
} // namespace corsika::process::observation_plane
diff --git a/Processes/ObservationPlane/testObservationPlane.cc b/Processes/ObservationPlane/testObservationPlane.cc
index 1df725253..957a06b97 100644
--- a/Processes/ObservationPlane/testObservationPlane.cc
+++ b/Processes/ObservationPlane/testObservationPlane.cc
@@ -41,7 +41,7 @@ TEST_CASE("ContinuousProcess interface", "[proccesses][observation_plane]") {
Vector<units::si::SpeedType::dimension_type> vec(rootCS, 0_m / second, 0_m / second,
-units::constants::c);
Line line(start, vec);
- Trajectory<Line> track(line, 10_m / units::constants::c);
+ Trajectory<Line> track(line, 12_m / units::constants::c);
// setup particle stack, and add primary particle
setup::Stack stack;
@@ -64,14 +64,14 @@ TEST_CASE("ContinuousProcess interface", "[proccesses][observation_plane]") {
Plane const obsPlane(Point(rootCS, {0_m, 0_m, 0_m}),
Vector<dimensionless_d>(rootCS, {0., 0., 1.}));
- ObservationPlane obs(obsPlane, "particles.dat");
+ ObservationPlane obs(obsPlane, "particles.dat", true);
obs.Init();
- [[maybe_unused]] const LengthType length = obs.MaxStepLength(particle, track);
- [[maybe_unused]] const process::EProcessReturn ret =
- obs.DoContinuous(particle, track);
+ const LengthType length = obs.MaxStepLength(particle, track);
+ const process::EProcessReturn ret = obs.DoContinuous(particle, track);
- SECTION("steplength") { REQUIRE(length == 10_m); }
+ REQUIRE(length / 10_m == Approx(1).margin(1e-4));
+ REQUIRE(ret == process::EProcessReturn::eParticleAbsorbed);
/*
SECTION("horizontal plane") {
@@ -81,16 +81,18 @@ TEST_CASE("ContinuousProcess interface", "[proccesses][observation_plane]") {
*/
}
- SECTION("inclined plane") {
+ SECTION("inclined plane") {}
+
+ SECTION("transparent plane") {
Plane const obsPlane(Point(rootCS, {0_m, 0_m, 0_m}),
- Vector<dimensionless_d>(rootCS, {1., 1., 0.5}));
- ObservationPlane obs(obsPlane, "particles.dat");
+ Vector<dimensionless_d>(rootCS, {0., 0., 1.}));
+ ObservationPlane obs(obsPlane, "particles.dat", false);
obs.Init();
- [[maybe_unused]] const LengthType length = obs.MaxStepLength(particle, track);
- [[maybe_unused]] const process::EProcessReturn ret =
- obs.DoContinuous(particle, track);
+ const LengthType length = obs.MaxStepLength(particle, track);
+ const process::EProcessReturn ret = obs.DoContinuous(particle, track);
- SECTION("steplength") { REQUIRE(length == 12_m); }
+ REQUIRE(length / 10_m == Approx(1).margin(1e-4));
+ REQUIRE(ret == process::EProcessReturn::eOk);
}
}
--
GitLab