diff --git a/corsika/detail/modules/ObservationPlane.inl b/corsika/detail/modules/ObservationPlane.inl
index 4229207f43de9004d88d920044ea1d92e8345dfe..42dd6ad5f4e2557dfe9ece6260d58f6da0a1ac06 100644
--- a/corsika/detail/modules/ObservationPlane.inl
+++ b/corsika/detail/modules/ObservationPlane.inl
@@ -16,49 +16,73 @@ namespace corsika::observation_plane {
std::string const& filename, bool deleteOnHit)
: plane_(obsPlane)
, outputStream_(filename)
- , deleteOnHit_(deleteOnHit) {
+ , deleteOnHit_(deleteOnHit)
+ , energy_ground_(0_GeV)
+ , count_ground_(0) {
outputStream_ << "#PDG code, energy / eV, distance to center / m" << std::endl;
}
- corsika::ProcessReturn ObservationPlane::doContinuous(
- corsika::setup::Stack::ParticleType const& particle,
- corsika::setup::Trajectory const& trajectory) {
+ ProcessReturn ObservationPlane::doContinuous(
+ corsika::setup::Stack::particle_type& particle,
+ corsika::setup::Trajectory& trajectory) {
TimeType const timeOfIntersection =
- (plane_.GetCenter() - trajectory.GetR0()).dot(plane_.GetNormal()) /
- trajectory.GetV0().dot(plane_.GetNormal());
+ (plane_.getCenter() - trajectory.getStartPoint()).dot(plane_.getNormal()) /
+ trajectory.getVelocity().dot(plane_.getNormal());
- if (timeOfIntersection < TimeType::zero()) { return corsika::ProcessReturn::Ok; }
+ if (timeOfIntersection < TimeType::zero()) { return ProcessReturn::Ok; }
- if (plane_.IsAbove(trajectory.GetR0()) == plane_.IsAbove(trajectory.GetPosition(1))) {
- return corsika::ProcessReturn::Ok;
+ if (plane_.isAbove(trajectory.getStartPoint()) ==
+ plane_.isAbove(trajectory.getPosition(1))) {
+ return ProcessReturn::Ok;
}
- outputStream_ << static_cast<int>(corsika::get_PDG(particle.GetPID())) << ' '
- << particle.GetEnergy() * (1 / 1_eV) << ' '
- << (trajectory.GetPosition(1) - plane_.GetCenter()).norm() / 1_m
+ const auto energy = particle.getEnergy();
+ outputStream_ << static_cast<int>(corsika::get_PDG(particle.getPID())) << ' '
+ << energy / 1_eV << ' '
+ << (trajectory.getPosition(1) - plane_.getCenter()).getNorm() / 1_m
<< std::endl;
if (deleteOnHit_) {
- return corsika::ProcessReturn::ParticleAbsorbed;
+ count_ground_++;
+ energy_ground_ += energy;
+ particle.erase();
+ return ProcessReturn::ParticleAbsorbed;
} else {
- return corsika::ProcessReturn::Ok;
+ return ProcessReturn::Ok;
}
}
- corsika::LengthType ObservationPlane::MaxStepLength(
- corsika::setup::Stack::ParticleType const&,
+ corsika::LengthType ObservationPlane::getMaxStepLength(
+ corsika::setup::Stack::particle_type const&,
corsika::setup::Trajectory const& trajectory) {
TimeType const timeOfIntersection =
- (plane_.GetCenter() - trajectory.GetR0()).dot(plane_.GetNormal()) /
- trajectory.GetV0().dot(plane_.GetNormal());
+ (plane_.getCenter() - trajectory.getStartPoint()).dot(plane_.getNormal()) /
+ trajectory.getVelocity().dot(plane_.getNormal());
if (timeOfIntersection < TimeType::zero()) {
return std::numeric_limits<double>::infinity() * 1_m;
}
- auto const pointOfIntersection = trajectory.GetPosition(timeOfIntersection);
- return (trajectory.GetR0() - pointOfIntersection).norm() * 1.0001;
+ auto const pointOfIntersection = trajectory.getPosition(timeOfIntersection);
+ auto dist = (trajectory.getStartPoint() - pointOfIntersection).getNorm() * 1.0001;
+ CORSIKA_LOG_TRACE("ObservationPlane::MaxStepLength l={} m", dist / 1_m);
+ return dist;
+ }
+
+ void ObservationPlane::showResults() const {
+ CORSIKA_LOG_INFO(
+ " ******************************\n"
+ " ObservationPlane: \n"
+ " energy in ground (GeV) : {}\n"
+ " no. of particles in ground : {}\n"
+ " ******************************",
+ energy_ground_ / 1_GeV, count_ground_);
+ }
+
+ void ObservationPlane::reset() {
+ energy_ground_ = 0_GeV;
+ count_ground_ = 0;
}
} // namespace corsika::observation_plane
diff --git a/corsika/modules/ObservationPlane.hpp b/corsika/modules/ObservationPlane.hpp
index d05e2b2762399c0b0acafe8fc1b0078ee2da894f..df57c77f4394f0fabecd44b0990d4cc828c8df4f 100644
--- a/corsika/modules/ObservationPlane.hpp
+++ b/corsika/modules/ObservationPlane.hpp
@@ -15,6 +15,7 @@
#include <corsika/setup/SetupTrajectory.hpp>
#include <fstream>
+#include <string>
namespace corsika::observation_plane {
@@ -29,17 +30,22 @@ namespace corsika::observation_plane {
ObservationPlane(corsika::Plane const&, std::string const&, bool = true);
void Init() {}
- corsika::ProcessReturn doContinuous(
- corsika::setup::Stack::ParticleType const& vParticle,
- corsika::setup::Trajectory const& vTrajectory);
+ corsika::ProcessReturn doContinuous(corsika::setup::Stack::particle_type& vParticle,
+ corsika::setup::Trajectory& vTrajectory);
- LengthType MaxStepLength(corsika::setup::Stack::ParticleType const&,
- corsika::setup::Trajectory const& vTrajectory);
+ LengthType getMaxStepLength(corsika::setup::Stack::particle_type const&,
+ corsika::setup::Trajectory const& vTrajectory);
+
+ void showResults() const;
+ void reset();
+ HEPEnergyType getEnergyGround() const { return energy_ground_; }
private:
- corsika::Plane const plane_;
+ Plane const plane_;
std::ofstream outputStream_;
bool const deleteOnHit_;
+ HEPEnergyType energy_ground_;
+ unsigned int count_ground_;
};
} // namespace corsika::observation_plane
diff --git a/tests/modules/CMakeLists.txt b/tests/modules/CMakeLists.txt
index 098412aaa8924debfa1bbbb92d1002835339dda8..335e68c7f456fde6e62ba416d70d0eabf1ce8f09 100644
--- a/tests/modules/CMakeLists.txt
+++ b/tests/modules/CMakeLists.txt
@@ -4,7 +4,7 @@ set (test_modules_sources
testStackInspector.cpp
#testTrackingLine.cpp
#testExecTime.cpp
- #testObservationPlane.cpp
+ testObservationPlane.cpp
testQGSJetII.cpp
# testSwitchProcess.cpp -> gone
# testPythia8.cpp
diff --git a/tests/modules/testObservationPlane.cpp b/tests/modules/testObservationPlane.cpp
index 63f5ea33768319f1fc7dc27a54f06177650d48a1..c1071a6ef324aa3fb65e5d0dad87304feb5163ee 100644
--- a/tests/modules/testObservationPlane.cpp
+++ b/tests/modules/testObservationPlane.cpp
@@ -22,7 +22,7 @@ using namespace corsika;
TEST_CASE("ContinuousProcess interface", "[proccesses][observation_plane]") {
- auto const& rootCS = RootCoordinateSystem::getInstance().GetRootCoordinateSystem();
+ auto const& rootCS = get_root_CoordinateSystem();
/*
Test with downward going 1_GeV neutrino, starting at 0,1_m,10m
@@ -38,19 +38,17 @@ TEST_CASE("ContinuousProcess interface", "[proccesses][observation_plane]") {
// setup particle stack, and add primary particle
setup::Stack stack;
- stack.Clear();
+ stack.clear();
{
auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
return sqrt((Elab - m) * (Elab + m));
};
- stack.AddParticle(
- std::tuple<Code, HEPEnergyType, corsika::MomentumVector, Point, TimeType>{
- Code::NuMu, 1_GeV,
- corsika::MomentumVector(rootCS,
- {0_GeV, 0_GeV, -elab2plab(1_GeV, NuMu::mass())}),
- Point(rootCS, {1_m, 1_m, 10_m}), 0_ns});
+ stack.addParticle(std::make_tuple(
+ Code::NuMu, 1_GeV,
+ corsika::MomentumVector(rootCS, {0_GeV, 0_GeV, -elab2plab(1_GeV, NuMu::mass)}),
+ Point(rootCS, {1_m, 1_m, 10_m}), 0_ns));
}
- auto particle = stack.GetNextParticle();
+ auto particle = stack.getNextParticle();
SECTION("horizontal plane") {
@@ -58,15 +56,15 @@ TEST_CASE("ContinuousProcess interface", "[proccesses][observation_plane]") {
Vector<dimensionless_d>(rootCS, {0., 0., 1.}));
ObservationPlane obs(obsPlane, "particles.dat", true);
- const LengthType length = obs.MaxStepLength(particle, track);
+ const LengthType length = obs.getMaxStepLength(particle, track);
const ProcessReturn ret = obs.doContinuous(particle, track);
- REQUIRE(length / 10_m == Approx(1).margin(1e-4));
- REQUIRE(ret == ProcessReturn::ParticleAbsorbed);
+ CHECK(length / 10_m == Approx(1).margin(1e-4));
+ CHECK(ret == ProcessReturn::ParticleAbsorbed);
/*
SECTION("horizontal plane") {
- REQUIRE(true); // todo: we have to check content of output file...
+ CHECK(true); // todo: we have to check content of output file...
}
*/
@@ -79,10 +77,10 @@ TEST_CASE("ContinuousProcess interface", "[proccesses][observation_plane]") {
Vector<dimensionless_d>(rootCS, {0., 0., 1.}));
ObservationPlane obs(obsPlane, "particles.dat", false);
- const LengthType length = obs.MaxStepLength(particle, track);
+ const LengthType length = obs.getMaxStepLength(particle, track);
const ProcessReturn ret = obs.doContinuous(particle, track);
- REQUIRE(length / 10_m == Approx(1).margin(1e-4));
- REQUIRE(ret == ProcessReturn::Ok);
+ CHECK(length / 10_m == Approx(1).margin(1e-4));
+ CHECK(ret == ProcessReturn::Ok);
}
}